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ABSTRACT 

We employ a flexible Bayesian technique to estimate the black hole mass and Eddington ratio 
functions for Type 1 (i.e., broad line) quasars from a uniformly-selected data set of ~ 58, 000 quasars 
from the SDSS DR7. We find that the SDSS becomes significantly incomplete at M B h < 3 x 10 8 M Q 
or L/LEdd ;$ 0.07, and that the number densities of Type 1 quasars continue to increase down to 
these limits. Both the mass and Eddington ratio functions show evidence of downsizing, with the most 
massive and highest Eddington ratio black holes experiencing Type 1 quasar phases first, although the 
Eddington ratio number densities are flat at z < 2. We estimate the maximum Eddington ratio of Type 
1 quasars in the observable Universe to be L/LEdd ~ 3. Consistent with our results in Paper I, we do 
not find statistical evidence for a so-called "sub-Eddington boundary" in the mass-luminosity plane of 
broad line quasars, and demonstrate that such an apparent boundary in the observed distribution can 
be caused by selection effect and errors in virial BH mass estimates. Based on the typical Eddington 
ratio in a given mass bin, we estimate typical growth times for the black holes in Type 1 quasars and 
find that they are typically comparable to or longer than the age of the universe, implying an earlier 
phase of accelerated (i.e., with higher Eddington ratios) and possibly obscured growth. The large 
masses probed by our sample imply that most of our black holes reside in what are locally early type 
galaxies, and we interpret our results within the context of models of self-regulated black hole growth. 

Subject headings: black hole physics — galaxies: active — quasars: general — surveys 



1. INTRODUCTION 

1.1. Background 

Understanding how and when supermassive black holes 
(SMBHs) grow is currently one of the great outstanding 
problems in extragalactic astronomy. The evolution of 
the galaxy and SMBH populations are not independent, 
as implied by established correlations between the mass 
of the SMBH, Mbh, and properties of the host galaxy 

"11991 
ver- 



sion (the Mbh~o~* relationshir 
120001 iMerritt fc Ferraresel l200lt 
concentration or Sersi c index 
Graham fe Driver] l2007f), bu lge 

19981: iMarconi fc Huntl l2003t lHaring fc Rixl 120041) . and 



. e.g.. iGebhardt et al.l 


Trcmaine et al. 
(Graham et al. 


2002), 


12001: 



gested that the origin of the scaling relationships emerges 
from the stochastic nature of the hierarchical assembly 
of black ho l e and stellar mass thro ugh galaxy mergers 
ljPenell2007t I.Tahnke fc Macci^l2lrrih . 

The galaxy and SMBH populations are also coupled 
because the fueling of SMBHs, which initiates AGN 
activity, depends on events that occur within and to the 
host galaxy. Many models have invoked major mergers 
of two gas-rich galaxies as the triggering mechanism for 
quasar activity that provide s the bulk of SMBH g rowth 
/e.g.. ISanders et all fl988t ISanders fc Mirabe| [1995 
' Kauffmann fc Haehneltl 120001: IWvithe fc Loebf booa 
Springe! et al.ll2005t iDi Matteo etall 2005; Si iacki et ail 
2007t iDi Matteo et al l 120081: pbpkins et all 120081: 



Somerville et all 120081: Tshenl I2009H Models have also 



bindin g energy (|Aller fc RichstoneT 20071 : IHopkins et all 
l2007a| ). Motivated by these empirical trends, a 
number of authors have invoked AGN 3 feedback as 
a means of regulating the SMBH's growth, which 
ties Mbh to properties of the host galaxy (e.g. , 
Silk & Recs 1998; Fabian 1999; Begel man fc N ath 2005; 
Murray et a l. 2005; Di Matt eo et al 112 005: Soringcl et al. 
20051: IHopkins et all l2006bt I Johansson et al.l l2009f ). In 
addition to regulating the growth of SMBHs, AGN feed- 
back has been invoked as a means of q uenching the 
growth of the most mas sive galaxies (e.g., iBower et all 
120061: ICroton et al.ll2006ll Alternatively, it has been sug- 
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invoked large-scale secular instabilities as a means of 
fueling AGN activity and growing SMBHs in disks (e.g.. 
Bower et al. 2006; Bournaud et~al]|2011llFanidakis et all 



20111) . At lower luminosities or redshifts, other fueling 
mechanisms may dominate the trigge ring of AGN 
activity (e.g., IHopkins fc Hernauistl |2009fl such as ex- 
terna l interactions (e.g. , ISerber et al.ll2006 ; lAlonso et all 
20071 IWjods fc Gelled 120071: iSilverman et al.l l20TlF 



Ellison et a. 



of gas (e.g. 



...20111: iLiu et al. 20121). sto chastic~accretion 
Hop kins fc Hernquistll2006l) . bar instabilities 



e.g.. IShlosman et all 119891: iGarcia-Burillo et al.l 120051: 
Hopkins fc Qua tacrt 2 0101), or stellar mass loss (e.g. , 
' Norman fc Scovilld 119881: ICiotti fc Ostrikerl Il997l 120071 : 
Kauffmann fc Heckmanl 120091: iHol I2009D . Regardless 
of the fueling mechanism, if the SMBH's growth is 
self-regulated via AGN feedback then the final mass of 
the SMBH after a fueling event is set by the binding 
energy of the bulge (jYounger et al.ll2008ft . 

While it is established that the galaxy and SMBHs 
populations are not completely independent, recent ob- 
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servational work has painted a more complicated picture 
regarding the connection between these two populations. 
In general, it is unclear if the scaling relationships be- 
tween the host galaxy and Mbh extend beyond classical 
bulges. T he scatter in the Mbh~o* relation is larger 
for spirals (|Grahamll200l iGultekin et al.ll2009fl . and ap- 
pears to increase at low Mbh such that most SMBHs li e 
below the Mbh~v* relation, (e.g.. iGreene et ai1l2010al ). 
Several authors have found evidence that there is no sig- 
nificant c orrelation between Mbh and pseudobulge lu- 
" Greene et all [20081: iGrahaml 1 2008: 



mmosity 



Graham fc Lil 



Jiang et al.l 120111: IKormendv et al.1 



20111) . Recently, IKormendv et al.l (|201lD argue that 
Mbh also does not correlate wit h properties of galaxy 
disks, while iGraham et al.1 (|2011[ ) concluded that Mbh 
does correlate with cr* for barred galaxies, even though 
their Mbh-&* relationship is offset from that of non- 
barred galaxies. These studies imply that the fueling 
mechansims of AGN activity, as well as the role of feed- 
back or other mechanisms for regulating the mass of the 
SMBH, are diverse and not necessarily uniform for all 
galaxy types. 

Ideally, in order to constrain models of SMBH fueling, 
growth, and impact on the host galaxy one would like 
to follow across time the stochastic process that is the 
evolution of SMBHs and their hosts. However, this is 
not possible, so demographic studies must be used to re- 
construct the evolution of these populations. Because of 
this, the demographics of SMBHs, AGN, and their host 
galaxies is one of the primary empirical tools that we have 
for placing constraints on astrophysical models of SMBH 
growth and fueling. Traditionally, studies of AGN de- 
mographics have focused on the luminosity function (e.g., 
for some recent estimates seelBovle et al.|2000l:lFan et a l.l 
200lHUeda et alll200l IWolf et al.l|2003HFan et afll20(> 
Barger et al. 2005; Hasingcr ct al. 2005; Richards et al 
200.1. 201 iFontanot et al| ^ (20071: iBongiorno et al 



2007 



Hopkins c t all 12007 



Croom et al.l 120091: I Jiang et al.1 120091: lAird et all 



Silverman et al.1 



2008 



2010; 



Willott et all I2010bb iFiore et all 120121: IShen fc Keliv 
2012 1 ) and their spatial c lustering (e.g.. iPorciani et al 
2004 [Croom et all 120051 : [Coil et all 120071: iMvers et al 
20071 IShen et al l l2007bl : Ishen et all 120091: IRoss et"a! 



20091: Ida Angela et al.l [20081 iMandelbaum et all l200l 



ng(._, 

Hickox et all 120091 120111: iWhite et all I2012D . From 
these studies, two important results have emerged. 
Many recent luminosity function studies lead to the 
discovery that the comoving number densities of 
more luminous AGN peak at earlier times than do 
the number den s ities of less luminous AGN (e.g. , 
Cowie et al.l 120031: ISteffen et all 12003b lUeda et all 120031: 



Hasinger et al.ll2005HBongiorno et al . 20071: iCroom eTall 
20091: iRigbv et al.ll2011[) . a phenomenon termed 'cosmic 



downsizing'. That the more luminous AGN population 
turns-off at earlier cosmic epochs implies that the more 
massive SMBHs grow first. The clustering studies have 
concluded that AGN at all redshifts typically live in dark 
matter halos of M h ~ 3 x 10 12 M Q , implying that AGN 
at higher z live in more biased environments. However, 
the halo mass can vary depending on t he AGN selection 
method (e.g.. Iffickox et alll2009l I2011D . 

More recently it has become possible to study 
the demographics of AGN host galaxy morpholo- 



gies out to z ~ 1-2, which provides insight into 
the fueling mechanisms of AGN activity. Stud- 
ies of optically-selected luminous quasars have found 
that they reside predominantly in elliptical galax- 
ies (e.g., [Kauffmann e jap 120031 : iDunlop et all 120031 : 
iKotilainen et alll2007l ). However, studies of X-ray se- 
lected AGN at lower luminosities (Lboi < a few x 
10 45 ergs _1 ) have found that AGN live in both elliptical 
and disk-dominated galaxies, and that their host galax- 
ies are no more likely to show evidence of disturbances 
or interaction signatures than their inactive counterparts 
(iGrogin et ^^O^^evce et alll2007l:IGabor et~aTll2009l: 



Cisternas et al 2011; 



Schawinski et allfeOllHPovic et al 



.20121; iKocevski et al. I I2012D . In contrast, iKoss et al 
(|2010l) find that hard X-ray selected AGN do pref- 
erentially reside in disturbed host galaxies at z < 
0.05. Studies of AGN host galaxy colors have 
found results consistent with the morphology stud- 
ies: namely, that AGN live in both early and late 
type galaxies (e.g.. iNandra et alll2007t [Coil et all 120091; 

schawinski et al.l 120091 ; IXue et al l 
In addition, lAird et all 



2009; 



Silverm an et al 
2010t iGeorgakakis et al.ll2011l 
(2012) find that the fraction of galaxies hosting an AGN 
is independent of stellar mass at fixed Eddington ratio. 
These results, particularly the large fraction of AGN in 
disk galaxies, imply that there is not one single dominant 
mechanism that initiates AGN activity and fuels SMBH 
growth, and that while major mergers are the most likely 
mechanism to fuel luminous quasars, they may not be the 
most frequent mechanism to fuel low luminosity AGN. In 
particular, the large fraction of AGN found in disk galax- 
ies implies that secular evolution plays an important role 
in fueling AGN activity. Unfortunately, there are dif- 
ficulties in identifying mergers as the trigger of AGN 
activity due to a possible time-lag betwe en the merger 
and the initiation of AGN ac tivity (e.g., iHopkins et all 
l2006al ISchawinski et al.ll2010D . combined with the rapid 
fading o f the morphologic al features indicative of a recent 
merger ()Lotz et alll2010D . 

Studies of AGN host galaxy demographics have 
also found evidence that star-formation rate and 
SMBH accretion rate, as traced through the AGN 
luminosity, are correlated both o n a cosmologica l 
level (e.g., ISilverman eFall 120081; lAird et all 120101) 
and within individual galaxies (e.g.. Jlmanishi &; Wada 
1 20041: (Schweitzer et alll2006l:IShi et alll2007l;lNetzer et al 
l2007bHLutz et al.ll2008tTShi et al.M2009t IMor et al.ll2012D . 
These results imply that SMBH fueling and growth 
is at least in part related to the presence of star 
formation, which is perhaps not surprising as both 
can be initiated by the presen c e of cold gas. Re- 
cently, iDiamond-Stanic fc Riekd (|2012f ) found that in 
local Seyfert galaxies the correlation between star for- 
mation rate and SMBH accretion rate is significantly 
stronger when limiting the analysis to the nuclear star 
formation rate, in agreement with results on SMBH fuel- 
ing obtai ned from smoothed particle hydrodynamic sim- 
ulations (Hop kins fc Quataertll2010l ). 

Another major advancement in AGN demographic 
studies has been the development of so-called 
virial mass estimates fo r Mbh (e.g., Wandc fet all 
19991: iMcLure fc Jarvis] 12001 iVestergaardl 12001 



Vestergaard fc Peterson! |2006[ ). This has opened up 



the possibility of studying AGN demographics with 



QUASAR DEMOGRAPHICS II 



3 



respect to Mbh- Compared to luminosity, studies of 
AGN demographics with respect to Mbh have the 
advantage that mass is a more stable quantity, in that 
Mbh can only increase and does so during active phases 
or through SMBH mergers. As a result, demographics of 
AGN M B h directly probe the growth of SMBHs. These 
virial mass estimates are derived by using the width of 
the broad emission lines as a proxy for the velocity dis- 
persion of the clouds emitting the broad emission lines, 
and the luminosity as a proxy for the broad line region 
size (iKaspi et al J 120051: IBentz et alll2009al: [Greene et al.l 
l2010bf K The virial mass estimates have a statistical 
scatter about the mass estimates d erived from rever- 
beratio n mapping dPeterson et al.l 12004 IBentz et all 
2009H) of ~ 0.4 d ex (e.g.. IVestergaard k PetersoiJl2006t 



Park et all |2012[ ). a lthough there may be add i tional 
systematic erro r s (iKrolikl |2001| ICollin et all |200€ 



Shen et all 120081: iFine et al.l 120081: [Marconi et all 1200 



Dennev et alJ 120091 : iRafiee k Halll l2011at ISteinhardd 



201 lh . The virial and reverberation mapping mass 
estim ates are calibrated to the local Mbh-v* relation- 
ship ([Onken et alj|2004t IWoo et all [20101: iGraham et all 

1ml). 

Numerous recent studies have used the virial 
mass estimates to study the demographics of AGN 
Mbh and Eddington ratio and their evolution (e.g.. 
Woo k Urrvll2002t : IVestergaardll200i iMcLure k Dunlcd 



200C 



2008 



2004 iKollmeier et all 1200a ISulentic et all 
Babic et al.l 120071 : iNetzer et alJ l2007at IShen et al.l 
Gavignaud et all 120081: IFine et all 120081: iTrump et al l 
20091 120111: iTrakhtenbrot et all I201U IRafiee k Halll 
2011b). In particular, the use of virial mass estimates has 



led several authors to estimate the black hole mass func- 
tion (BHMF) and Ed dington ratio func t ion (BHERF) 
for T y pe 1 quasars, dWang et all 12001 iGreene k Hoi 
20071: IVestergaard et al. 2008; Vcstcrgaar d k Osmerl 



20091 IKellv et all 120091 iSchulze k Wiso tekl l2010l 



Kelly et all 1201a iWillott et all I20l6at IShen k KeTivl 
2010,12013), which gives the comoving number density of 



Type 1 quasars as a function of Mbh an d L/Lsdd- The 
BHMF and BHERF are important because they provide 
a complete census of the active SMBH population with 
respect to Mbh and L/LEdd-, making them one of 
the primary empirical tools for constraining models 
of SMBH growth. Indeed, many models for SMBH 
growth have made predictions for t he BHMF and 
BHERF at a variety of red s hifts (e.g., [C attanc o et all 
20051 : iDi Matteo et alj|200l [Hopkins et al.l [2008; Slier 



Tanaka k Haimanl 120091: IVolonteri k Begelmar 
20iq |Fanidakis et al. 2011, 20121: iNataraian fc Volonteril 
2012t iDraper k Ballantvnd 120121 ). In addition, M B h 
is a fundamental physical parameter of black hole 
accretion flows, making the BHMF and BHERF im- 
portant for studies of accretion physics, as it describes 
which regions of the Mbh~L/ LEdd plane are probed 
by current and future surveys. Based on AGN num- 
ber densities, several groups have found evidence for 
down sizing in Mbh of SMBHs in Type 1 quasa rs 
(e.g.. IVestergaard k Osmerl l2009l lLabita et al.l l2009allbl: 
IKellv et alj|2010t IShen k Kellvll2012D . implying that at 
least some of the downsizing in the AGN luminosity 
function may be driven by downsizing in Mbh- 

Early estimates of the BHMF were obtained by directly 
binning up the virial mass estimates and applying the 



}/V m ax correction (TWang et alj2006tlG reene "&Hol2007l: 
IVestergaard et all 120081: IVestergaard k Osmerl 120091) , 
where the 1/V max correction is the same as that used 
to estimate the luminosity function. However, BHMFs 
obtained in this manner suffer from both incompleteness 
due to the sample flux limit and artificial broadening 
due to the statistical error in the virial mass estimates 
(IKellv k Bechtoldl [20071: IShen et all [20081 IKellv et all 
2009). The incompleteness arises because there is a 
large range in luminosity at fixed Mbh, scattering some 
quasars at fixed Mbh above the flux limit and some be- 
low. The artificial broadening arises because the sta- 
tistical error in the virial mass estimates scatters more 
quasars into bins of higher Mbh than lower Mbh when 
the BHMF declines toward high er values of Mbh - In 
order to account for these effects, IKellv et al.l ([20091 ) de- 
veloped a Bayesian technique for estimating both the 
BHMF and BHERF of Type 1 quasars which corrects 
for incompleteness and the effects of statistical error in 
the virial estimates; they used their method to estimate 
the lo cal BHMF of Type 1 quasars. Sch ulze k Wisotzkl 
(2010) developed and used a similar method to estimate 
the local BHMF and BHERF of Type 1 quasars, al- 
though they did not correct for th e error in the virial 
mass estimates. IKellv et al.l (I2010L here after K10) ex- 
tended the method of IKellv et all " ?|2009h to use a more 
flexible Eddington ratio distribution at fixed Mbh to es- 
timate the BHMF at z > 1 from the virial mas s estimates 
from the SDSS DR3 (IVestergaard et al.ll2008l) . K10 con- 
cluded that there is a broad range in L/LEdd for Type 1 
quasars, that most Type 1 quasars are not radiating near 
the Eddington limit, and that their results are consistent 
with models for self-regulated black hole growth. 

To study SMBH growth, one would ideally like to 
estimate the BHMF and BHERF for all SMBHs, i.e., 
including both active (broad line or not) and inac- 
tive SMBHs. Unfortunately, this is currently only 
possible beyond the local universe using techniques 
which indirectly estimate the BHMF. Several authors 
have derived local BHMFs for all SMBHs from indi- 
vidual sources using the scaling relations hips between 
Mbh and host spheroidal pr o perties (e.g.. ISalucci et all 
2Q0j IWfcjTremaind l200l lA llcr k Richst ond 120021 
Marconi et all 12004 iShankar etaLI 12004 [ Tundo et al 
lYu k Lul |2008rlShankar et all 120091 IVika et al 



2007 



20091 ). Many studies have found evidenc e for evo- 
lution in the scaling relationships (e.g., iTreu et all 



200llPeng et al.ll200a iTreu et all[2007l lWoo et al.H200 



Decarli et al.l 120101: iMerloni et all 120101: iBennert et all 
20101 but see Lauer et al. 2007, Shen k Kelly 2010, 



Schulze k Wisotzki 2011, and Portinari et al. 2012 for 
cautionary notes), and, while there have been attempts 
to estimate the BHMF based on an ass umed form of 
the evolution in the scalin g relationships (jTamura et all 
12001 iShanker et all [20091 iLi et al.l [20T1 . the quanti- 
tative form of the evolution has not been sufficiently 
precise to motivate widespread use of the scaling re- 
lationships to estimate BHMFs beyond the local uni- 
verse. Instead, indir ect techniqu es based on variations 
of the argument of (|Soltanl [l982f ) are used to estimate 
the BHMF for all SMBHs beyond the local universe. 
The basic idea behind these techniques is to step back- 
ward in time from the local BHMF using the AGN lu- 
minosity function as a probe of how fast the BHMF 
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is changing in time. A number of groups have used 
variations on this te chnique (e.g., lYu k Tremaind 120021 : 
Marconi et al.112004 lMerlonill2004[Hopkins et al.ll2007bl: 
Merloni k Heind 120081: IShankar et all 120091: ICaol l2010h . 
and have generally concluded that most, if not all, of the 
local BHMF can be explained as the relic of AGN ac- 
tivity, with SMBH growth being dominated by periods 
when the SMBH was radiating near the Eddington limit. 
In addition, these studies have generally concluded that 
SMBH growth is anti-hierarchical, in agreement with 
the cosmic downsizing seen in the AGN luminosity func- 
tion studies. For a rev iew of the BHMF of SMBHs, see 
iKellv k Merloni ([2011 1. 

It is currently only possible to obtain estimates of Mbh 
for large samples of individual sources out to z > 1 us- 
ing the virial mass estimates. Therefore, while Type 
1 quasars represent only a fraction of the SMBH pop- 
ulation, they are the only objects for which it is cur- 
rently possible to obtain an estimate of the BHMF be- 
yond the local universe that is derived using estimates 
of Mbh from individual sources. Moreover, the BHMF 
and BHERF of Type 1 quasars provides the distribution 
of Mbh and L/L,Edd for a subset of the population of 
SMBHs that are actively growing at any given redshift. 
The discovery of Type 1 quasars with Mbh ~ 10 9 A^ffr , 
out to z ~ 6-7 Oiling et alJ l2007t iKurk et al.1 120071 : 
IWillott et aTl feOlOa; M ortlock et all 12011( 1 implies that 
the BHMF and BHERF of high-redshift Type 1 quasars 
places strong constraints on models for the formation and 
growth of SMBH seeds (e.g.. lHaiman k Loebl[200l . 

1.2. Comparison with \Shen & KeM 1201& ) 

This is the second paper in a two-part series to study 
the demographics of Type 1 quasars out to z ~ 4.75 in 
the two-dimensional quasar mass and luminosity space. 
Such studies represent a next step in Type 1 quasar de- 
mographic studies by placing joint constraints on the dis- 
tributions of their instantaneous activity (i.e., luminosity 
function) and assembly history (i.e., BH mass function) 
via the Eddington ratio distribution. As such, the 2D 
demographic studies provide an advanced view of the 
cosmic growth of the SMBH population compared to us- 
ing the luminosity function or BHMF alone. Our study 
is motivated by both the recent availability of ~ 58, 000 
uniformly selected Type 1 quasar s from the S PSS DR7 
with virial mass estimates from IShen et ail ((201 lh . as 
well as the recent advancements in statistical techniques 
for studying AGN demographics. The incredible size of 
this data set combined with our new powerful statistical 
tools provides us with an unprecedented ability to char- 
acterize the demographics of Type 1 quasars at a vari- 
ety of redshifts. The broad redshift range of our sample 
(0.4 < z < 4.75) allows us to probe the evolution of a 
subset of the actively growing SMBH population during 
the e pochs over whi ch ~ 95% of SMBH growth occured 
(e.g.. IShankar et al.ll2009l) . 

In the first paper of this series, IShen k Kellvl (120121 
hereafter Paper I), we presented our sample as well as 
the binned estimates of the luminosity function and the 
BHMF. We also extended the statistical method of K10 
to estimate the BHMF independently in different redshift 
bins, and to include a luminosity-dependent bias term for 
the virial mass estimates. The former improvement en- 
ables us to study the evolution of the BHMF and BHERF 



without assuming a parameteric form for this evolution, 
while the latter improvement enables us to study the ef- 
fects of a systematic error in the virial mass estimates 
with luminosity. A luminosity-dependent bias term is 
motivated by the small scatter in the virial mass esti- 
mates at fixed Mbh and L, as obs erved by several recent 
studies of AGN demograph i cs (iKollmeier et al.1 120061: 
I Shen et all 120081: iFine et all 120081: ISteinhardt k Elvisl 
I2010H IShen k Kellvl 120 10L K10L The results from our 
model implied that virial mass estimates derived from 
the FWHM of the Mgn and Civ emission lines system- 
atically overestimate Mbh at higher than average lumi- 
nosity, assuming that the statistical scatter in the mass 
estimates at fixed Mbh and L is constant. 

In Paper I we also used our model to study the luminos- 
ity function of Type 1 quasars, and presented the BHMF 
and BHERF derived from our model. We studied the 
luminosity function below the flux limit, assuming our 
model for the SMBH mass and Eddington ratio distri- 
butions. We found evidence for downsizing in both the 
luminosity function and BHMF. In addition, we found 
that our model luminosity function makes reasonable 
predictions when extrapolated to ~ 3 mag fainter than 
our SDSS DR7 sample, and the extrapolations were of- 
ten consistent with luminosity functions estimated from 
deeper surveys. We also found evidence that the aver- 
age Eddington ratio of Type 1 quasars increases towards 
higher redshift, in qualitative agreement with predictions 
from recent h ydrodynamic simula tions of SMHB growth 
at z > 4 (e.g.. lDeGraf et alJl20H l. 

In this paper we build upon the work of Paper I and 
focus our analysis and discussion on the BHMF, the 
BHERF, and quantities that can be derived from them. 
The main differences between Paper I and this paper are: 

• We make several improvements to our statistical 
model. We extend the statistical model of K10 
and Paper I to model the joint distribution of Ed- 
dington ratio and Mbh at a given redshift as a 
mixture of 2-dimensional log-normal distributions. 
This provides us more flexibility compared to Pa- 
per I, where we assumed that the Eddington ratio 
distribution at fixed z and Mbh could be described 
as a single log-normal distribution with geometric 
mean depending linearly on log Mbh- In addition, 
we also use a Student's i-distribution to model the 
distribution of measurement errors in the mass esti- 
mates, which downweights outliers compared to the 
previously-used Gaussian distribution. This latter 
improvement was incorporated to make our results 
more robust against outliers in the mass-luminosity 
plane, as such sources may be subject to uniden- 
tified systematic error. We made these improve- 
ments to test the robustness of the key conclusions 
from Paper I. 

• Because we have employed a more flexible model 
for the BHERF, the focus of this paper is on the 
BHMF and BHERF of Type 1 quasars. As such, 
in this paper we discuss the BHMF and BHERF in 
greater depth than we did in Paper I, and discuss 
the distribution of Eddington ratios at fixed Mbh- 
In addition, we discuss what the derived BHMF 
and BHERF imply with regard to the growth of 
supermassive black holes. 



QUASAR DEMOGRAPHICS II 



5 



• We focus our presentation and discussion on the 
BHMF and BHERF based on the usual assump- 
tion that the mass estimates are unbiased. This is 
the assumption that is typically made in the litera- 
ture, but is different from Paper I where we allowed 
the mass estimates to have a luminosity-dependent 
bias. In this paper we present the results for the 
model assuming the estimates are unbiased, but 
also discuss how allowing a luminosity-dependent 
bias changes the results. Between Paper I and this 
paper we cover both simple models for the behav- 
ior of the virial mass estimates with respect to the 
BHMF and BHERF. 

• Unlike in Paper I, we do not present or discuss the 
optical luminosity function implied by our model, 
as there is little difference from Paper I. Instead, 
we focus on the BHMF and BHERF. 

• Paper I goes into greater depth with respect to bi- 
ases in the virial mass estimates, and what we can 
conclude about such biases based on AGN demo- 
graphics. That is not the focus of this paper. 

Between Paper I and this paper we present an in-depth 
analysis of Type 1 quasar demographics with respect to 
the properties of their virial mass estimates, their lu- 
minosity function, their black hole mass function, their 
black hole Eddington ratio function, and their behavior 
in the mass-luminosity plane. 

This paper is organized as follows. We describe our 
statistical model in §S1 We present our estimated BHMF 
in S}3] and BHERF in and discuss our results in <|5j 
We summarize our results in Ej6j Throughout the pa- 
per we adopt a flat ACDM cosmology with cosmologi- 
cal parameters ft\ — 0.7, ilo — 0.3, h = 0.7, to match 
most of the recent quasar demographics studies. Vol- 
ume is in comoving units unless otherwise stated. We 
distinguish virial masses from true masses with a sub- 
script v ; r . For simplicity, quasar luminosity is expressed 
in terms of the rest-frame 2500 A continuum luminosity 
(L = XL\). The conversion between the 2500 A con- 
tinuum luminosity and the absolute i-band magnitude 
n ormalized at z = 2 (M j\z = 2]) is given by eqn. (4) 
in Richar ds et al.l ((2006). Lower case letters refer to log- 
arithms of quantities based on mass or luminosity, e.g., 
I = logi, m B H = log M B h- 

2. THE STATISTICAL MODEL AND POSTERIOR 
DISTRIBUTION 

In this work we expand on the statistical model de- 
scribed in Paper I, which its elf is an expansio n of the 
statist ical model d e velope d bv lKellv et al.l ()2009l . see also 
K10). IKellv et al.l (|2009D developed a Bayesian method 
for estimating the BHMF and Eddington ratio distri- 
bution from the broad line mass estimates that cor- 
rects for both the statistical scatter in the mass esti- 
mates and incompleteness. They modeled the distribu- 
tion of Type 1 quasars in the mass-redshift plane as a 
mixture of log-normal distributions, and the Eddington 
ratio distribution at fixed Mbh as a single log-normal 
distribution whose geometric mean depended linearly on 
log Mbh- They assumed tha t the mass estimates were 
unbiased. IKellv et ail (|2010l) expanded this model and 
used a mixture of log-normals for the Eddington ratio 



dist ribution at f i xed M bh- In both IKellv et al.l (|2009l ) 
and IKellv et all (|2010| ) the Eddington ratio distribution 
was assumed to not evolve. In this section we describe 
our expansion to the model of Paper I and summarize 
the important aspects of our statis tical model; furth er 
details can be found in Paper I and Kell y et all (J2009) . 

2.1. Mixture of Log-Normal Functions Model 

Motivated by the observed small statistical scatter in 
the mass esti mates for the SPS S, Paper I expanded on 
the model of IKellv et all (2009) to incorporate a more 
flexible model for the error distribution of the mass esti- 
mates. In both Paper I and this paper the error distribu- 
tion is modeled as a Gaussian distribution with unknown 
variance and optionally and unknown mean. Because 
we use FWHM-based virial mass estimates, our model is 
only with respect to these mass estimates; mass estimates 
based on the line dispersion or other line width measures 
may have a different error distribution. Paper I assumed 
that the mass estimates were unbiased at the average 
luminosity as a function of Mbh, but that the mass esti- 
mates potentially exhibited a luminosity-dependent bias 
at fixed mass. They assumed the following model for the 
mass estimates: 

m mr = m B H + P[l - E(l\m B H)} + o m ie m \. (1) 

Here, m vir = log M vir ,m B H = log M B h, I = logi, 
E(1\itibh) is the expectation value of log L at fixed Mbh 
and e m i is a random variable drawn from the standard 
normal distribution. The term f3 models how the sys- 
tematic error (i.e., the luminosity-dependent bias) in the 
mass estimates scales with luminosity, and a m i is the 
standard deviation in the mass estimates at fixed Mbh 
and L. The standard mass estimates assume /3 = and 
a m i ~ 0.4 dex. 

The motivation for modeling the mass estimates ac- 
cording to Equation (TT]) is illustrated in Figure [TJ Here 
we show the distribution of the error in the mass es- 
timates at fixed true mass as a function of luminosity, 
assuming a value of (3 = 0.5 and a m i = 0.25. When av- 
eraging over a broad range of luminosity, as is done for 
the reverberation mapping sources, the distribution of 
the mass estimate errors is broad. However, when limit- 
ing ourselves to a narrow luminosity range at the bright 
end, as in a flux-limited sample like the SDSS, the mass 
estimates exhibit a bias and smaller scatter. This model 
is thus one way of reconciling the fact that the scatter in 
the mass estimates for the SDSS imply smaller uncertain- 
ties in M v i r than does the scatter in the mass estimates 
for the reverberation mapping sample. Some physical 
reasons to expect such a luminosity-dependent bias were 
further discussed in §3.2.1 of Paper I. Alternatively, an- 
other possibility is that the mass estimates at fixed mass 
and luminosity are always unbiased (i.e., /3 = 0), but o~ m i 
decreases toward higher luminosity. This possibility is 
also illustrated in Figure [1] Currently it is not possible 
to unambiguously distinguish between these two possi- 
bilities solely from AGN demographics, and both situa- 
tions may be at work. As discussed in Paper I, there is 
some indication from the reverberation mapping data for 
NGC 5548 that /3 > for FWHM-based virial mass esti- 
mates, although this result is only moderately significant 
at 2.4cr. As such, in this work we obtain constraints on 
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the BHMF and Eddington ratio distribution under each 
of these two models for the error in the mass estimates. 

The values of /3 and a m i are assumed to be different 
when different emission lines are used to calculate M V i r . 
In this work we use the H/3, Mgn, and Civ emission lines. 
These lines have different advantages and disadvantages, 
and the behavior of the mass estimates derived from the 
Mgn and Civ lines are less well understood than those 
derived from H/3. We refer the reader to Paper I for a 
discussion of the issues and possible biases surrounding 
virial ma ss estimates derive d from different emission lines 
(see also lShen fc Liull2012L and the discussion in § 15. 3|) . 

In this work we model the joint distribution of Mbh 
and L for Type I quasars in a given redshift bin as a 
mixture of K 2-dimensional log-normal distributions: 

p(m BH ,l\Tr,n,i:) = 

t ^17* e Xp {-i(x- Mfc fE-(x -, k )}. (2) 

Here, x = [tubm, 1} T ', x T denotes the transpose of x, tt — 
(tti, . . . ,tt k ), /u = {hx,...,hk), and £ = (Si, . . . , Y, K ). 
As in Paper I we model the distribution of z in each 
redshift bin as a Pareto distribution (i.e., a power-law): 



p(z\j) = 



(1+7)* 



1+7 
^max 



1+7 ' 



(3) 



In this work we use a value of if = 5, which provides 
considerable flexibility. There was little difference in em- 
ploying larger values of K and such larger values did not 
justify the increased computational burden. Note that 
the comoving number density per 3-dimensional box de- 
fined by (rriBH, h z) and (wish, I, z) + (dmsH, dl, dz) is 

(f>(m B H, l,z) = N ( ^— ) p(m B H,l\n, S)p(z|7), 



(4) 

where N is the total number of Type I quasars in the 
observable universe, and dV/dz is the derivative of co- 
moving volume with respect to redshift. The mass or 
luminosity function is obtained by integrating Equation 
l[4|) over I or ttibh, respectively. 

In the absence of measurement error, Equations ([!])- 
^ provide us with all of the necessary ingredients to 
construct the likelihood function; i.e., these equations 
tell us how to connect the mass and luminosity functions 
to the observed distribution of virial mass estimates, lu- 
minosities, and redshifts. However, in reality the line 
widths and luminosities are measured with error, and 
therefore so are the virial mass estimates. In general, 
the measurement error on the FWHM dominates over 
the measurement error on the luminosity, so we only 
account for measurement errors on the virial mass es- 
timates. Denote the measured virial mass estimate as 
M v i r (rh v ir = log M v i r ). In this work we model the mea- 
surement error distribution of the virial mass estimates 
using a Student's ^distribution centered at the true value 
of M v i r : 

P\Jfl"vir y^virj 

r((* + i)/2) . i / 

^mcasl 



1 + i 

v 



>+l)/2 



Here, v is the degrees of freedom of the ^distribution, 
fmcas is the fixed measurement error amplitude, which 
is calculated from the line fitting procedure, and T(-) is 
the Gamma function. We use the student's t-distribution 
because it is considered a robust alternative to the nor- 
mal distribution. In the limit v — > oo the t-distribution 
converges to the normal distribution, but for finite v 
the ^-distribution has heavier tails. As a result, the t- 
distribution downweights the influence of outlying values 
of the virial mass estimates, which may be caused by sys- 
tematics involving bad line width measurements, or due 
to the presence of a population of AGN for which the 
virial mass estimates are subject to a large systematic 
error. In this work we use v = 8, which is a typical value 
for robust analysis; there is little difference in using sim- 
ilar values of v. This is a slight improvement over Paper 
I, where we modeled the measurement error distribution 
as a normal distribution. 

For a given redshift bin, denote the set of parameters 
for our statistical model as 9 = (tt, fx, S,7, 0,cr m i), the 
vector of measured virial masses as M„i r , the vector of 
measured luminosities as L, the vector of measured red- 
shifts as z, and the n umber of Type 1 quasars in our 
sample as n. Following iKellv et all ((20091) , we derive the 
posterior distribution from Equations (HJ-© and ([S]) as 

p(6\M mr ,L,z)=p(9)[p(I = l\6)}- n 

n rOO 

x I1p( z j|7) / 

i)p(m V i rt i,k\d)dm v i r! i, 
(6) 

where 

p(m v ir,i,li\6) = 

K n k f I 

27r|Vfc| 1 /2 exp \ 2^™^ ~ Mfc) T ^ _1 (x«ir,i ~ Mfc) 



v k 



Var(m V i r \k) Cov(m V i r ,l\k) 
Cov(m V ir,l\k) T,L,k 



(7) 
(8) 

(9) 



(5) 



Var(m vtr \k) = T, M ,k + ft 2 {^L.k - kl^M,k) + cr^ 

(10) 

Cov{m vir ,l\k) = E M L,k + ^( s i,fe- s ifi,fc/ S JW,fe)- (11) 

Here, p{9) is the prior on 9, n is the total number of 
data points in the redshift bin, p(I = 1\9) is the proba- 
bility of a Type 1 quasar from the redshift bin of interest 
making it into our sample, Vk is the covariance matrix of 
log M V i r and logL for the fc th log- normal function, and 
^M,k,^L,k, and T,ML,k denote the variance in log-Me^, 
variance in logL, and covariance between log Mbh and 
logi for the k th log- normal function, respectively. The 
parameters, 9, are estimated independently in each red- 
shift bin. 

As in Paper I, we define L to be the luminosity 
at 2500A, which we derive from the ?-b and magnitude 
accord ing to the prescription given in Richard s et al.l 
(2006). The term p{I = 1|0) is calculated from the SDSS 
selection function, s(L, z), as 

/•oo />oo 

p(I=l\9)= / s(L, z) P {L\tt, h, E)p(z\i) dL dz. 
Jo Jo 

(12) 
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Fig. 1. — Illustration of the two models for the distribution of virial mass estimate errors used in this work. Mock data are generated using 
Eqn. |[T}. in both panels £?(log L\ Mgjj) denotes the mean value of logL at fixed Mgjj. The left panel shows a model with a luminosity 
dependent bias (/3 > 0) and constant cr m ; with luminosity, while the right panel shows a model with = and a decrease in the amplitude 
of the error in the mass estimates toward higher L. In both models the dispersion in virial estimates decreases when limited to only those 
quasars in the bright tail of the luminosity distribution. However, the virial mass estimates are also biased if the small dispersion in mass 
estimates observed in the SDSS is caused by the situation depicted in the left panel (/? > 0). 



The selection func tion of our sample is t he same as that 
for the sample in iRichards et al.l ((20 06) . given in that 
paper. 

2.2. The Prior Distribution 

Flexible models such as the mixture of log-normal func- 
tions model we use enable modeling of a broad range of 
distributions, but they can also suffer from overhtting 
because the data do not provide enough information on 
the structure of the distributions. It is often useful to 
impose priors invoking smoothness on the estimated dis- 
tributions, otherwise highly 'wiggly' BHMFs are consid- 
ered just as likely as smooth ones a priori. In addition, 
because our sample is truncated, there is little to no infor- 
mation from the data on the distribution in the Mbh~L 
plane below the flux limit. Therefore, it is necessary to 
impose several constraints on 9 via the prior distribu- 
tion in order to keep the solution from wandering into 
unreasonable regions. 

The prior constraints that we impose are as follows: 

• The standard deviation of each log- normal function 
must be between 0.2 and 1.0 dex for both log Mbh 
and log L/L Edd . 

• The mean of log Mbh for each log-normal func- 
tion must be between 6.0 and 10.0. This reflects 
our assumption that the BHMF must decrease at 
M BH < 10 6 M Q and M BH > 10 10 M Q . 

• The mean value of log L / LEdd for each log- normal 
function must be between -3.0 and 0.0. The 
lower bound was chosen because accretion flows 
are thought to undergo a state-transition near 
L/LEdd ~ 10 _2 -10 -3 , and the broad line region 
is not expected to exist below this c ritical Edding- 
ton ratio (e.g.. ICzernv et alj |200l iHopkins et all 
120091 : iTrump et al.l 1201 ll ). The upper bound was 
chosen to reflect our assumption that the BHERF 
must decrease above the Eddington limit. 

• The fraction of Type 1 quasars in each log-normal 
function radiating at log Lj L B dd < —3.5 must 
be less than 1%. This constraint was added in 



addition to the constraint on the mean value of 
log L I LEdd to ensure that the estimated distri- 
bution of Type 1 quasars declined strongly at 
/• /./:././:,. 10 ''. 

In order to test how sensitive our results are to these 
prior constraints, we have also performed the analysis 
allowing the standard deviation for both log Mb # and 
log L/LEdd to be as high as 2.0 dex, and extending the 
lower limit for the mean value of log L / LEdd- Extending 
the upper bound on does not significantly change 

our results beyond slightly increasing the uncertainties, 
with the only exception being that derived value for the 
maximum value of Mbh m the observable universe for 
a Type 1 quasar is very sensitive to this upper bound; 
this is discussed further in Section 13.31 In addition, the 
shape of the BHMF and BHERF at M BH > 3 x 10 8 M Q 
and L/LEdd ^ 0.05 are not sensitive to the lower bound 
on the mean value of log L / LEdd- However, the normal- 
ization of the BHMF and BHERF arc sensitive to this 
lower bound, as this lower bound is directly related to 
how many Type 1 quasars are sample is missing. Luck- 
ily this is not a problem for our analysis as none of our 
scientific conclusions depend on the absolue normaliza- 
tion of the BHMF and BHERF. 

Beyond the bounds listed above, we use a prior for the 
log-normal para meters that is based o n the 'partially- 
proper prior' of iRoeder fc W asscrman (1997). We con- 
strain the K log-normal functions by ordering their val- 
ues of (j,k by increasing mean luminosity, i.e., the k = 1 
log-normal function has the faintest mean luminosity, 
and the k = K log- normal function has the brightest. 
In addition, denote [ih.k and (iM,k to be the components 
of fik corresponding to the mean logL and log Mbh i re- 
spectively. Our prior for the luminosity components of 9 
is 

K 
fe=2 

(13) 

pfo.i) = ~* min fl- f £ - 1 - 1 (14) 

^min 'max \ *min 'max / 
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2Trh(k,k- 1) 



h L (k,k-l) 



x N(n Lyk \n L<k _x,hL{k, k - l))H(jiL,k - A*i,fe-i) (15) 



h L (k,k-l) 



J L,fc-l 



(16) 



Here, l m i n and / max are the minimum and maximum pos- 
sible values of fJ,L,k, h,L,(k,k — 1) denotes the harmonic 
mean of £z,,fc and Y>L,k-i, $(•) is the standard normal 
cumulative distribution function, N(x\fi, V) denotes a 
Gaussian function in x with mean /i and variance V, 
and H (•) is the Heaviside step function. Given our prior 
constraints, i min is the value of logL for M BH = 10 6 Af Q 
and L/LEdd = 10 -3 , and l max is the value of log!/ for 
Mbh = 10 10 Afo and L/LEdd = 1. We chose the prior on 
fiL,i because it is the probability distribution for the min- 
imum value of a sample of K random variables uniformly 
distributed between l m - ln and Z max - The prior on /i^.fe for 
k > 1 was chosen to enforce smoothness in the estimated 
luminosity function, as it places higher probability on 
solutions where the individual log-normal functions are 
close together with respect to their average standard de- 
viations. Our prior on was therefore constructed to 
not favor a particular value for the centroid of the lumi- 
nosity function, under the constraints listed above, but 
to favor smooth luminosity functions. 

We also place a similar 'smoothness' prior on /xm, but 
without enforcing a particular ordering of the log-normal 
functions: 

t | V -n TT 1 f 1 (hiM,k-m) 2 

p{p M \^M,m) = [[ /n i exp 



fc=i 



flM 



(17) 
(18) 



The 'hyper-parameter' m is the mean value of the {fiM,k} 
and is an additional parameter. We assume a uniform 
prior on rh over all possible values; note that even though 
we do not place explicit bounds on m, such bounds are 
implied by the bounds placed on jJ,M,k- We place uni- 
form priors on E^/. and Eju,t under the constraints given 
above. We obtain the prior on TiML,k by placing a uni- 
form prior on the angle that the slope of the mean value 
of log L/LEdd as a function of Mbh makes with the hor- 
izontal. For 7r, we assume a Dirichlet prior with param- 
eter a = 0.1. Finally, we use the same prior on (3 and 
<7 m ; as that described in § 3.2.1 of Paper I, which is con- 
structed to give results consistent with the reverberation 
mapping sample. 

3. THE BLACK HOLE MASS FUNCTION 

We applied our Bayesian method to our sample derived 
from the SDSS DR7 data set for both the model with a 
luminosity-dependent bias (/3 ^ 0) and the model with 
decreasing scatter in the mass estimates at higher lumi- 
nosity (/? = 0). For the latter model (/? = 0) we do not 
explicitly fit for a dependence of o m i on L, but rather 
just assume a single value of a m i over the luminosity 
range in each redshift bin. This should be an adequate 
approximation as the luminosity range of our data set 



in each redshift bin is narrow. Details of the sample are 
described in Paper I 

Assuming the model with a luminosity-dependent bias, 
we estimate the slopes of the bias to be /3np ~ 0.1 ± 
0.1,/3 Mg n ~ 0.39 ±0.13, and /3 C iv ~ 0.40 ± 0.11 for 
the H/3,Mgn, and Civ emission lines, respectively. For 
the f3 model we estimate the scatter in the mass 
estimates at fixed luminosity and BH mass to be uh/3 ~ 
0.24±0.03, o-Mgii ~ 0.25±0.02, and cr C i V « 0.21±0.01 for 
the H/3,Mgn, and Civ emission lines, respectively. These 
values arc significantly different from the traditionally 
assumed values of /3 = and <jbl — 0.4 dex. However, 
we note that at z > 3.5 the derived values for the Civ 
line become /3civ ~ 0.3 ± 0.3 and ctqiv ~ 0.39 ± 0.04, 
probably reflecting the broader range in luminosity at 
these redshifts due to the deeper SDSS flux limit. This 
increase in <j m i as the range in luminosity is increased 
supports our hypothesis that the dispersion in the virial 
mass estimate error is luminosity-dependent. 

These results on the luminosity-dependent bias model 
are similar as that obtained in Paper I using a simpler 
model for the Eddington ratio distribution, showing that 
many of the results of Paper I are robust against the as- 
sumed Eddington ratio distribution. We therefore con- 
clude from this that either the statistical scatter in the 
FWHM-based broad line mass estimates is correlated 
with luminosity, producing a luminosity-dependent bias, 
or that the scatter in the mass estimates is smaller for the 
luminous quasars probed by the SDSS, or some combina- 
tion of these two possibilities. In Paper I we have already 
discussed the implications of these results for the broad 
line mass estimates, and, as the focus of this paper is on 
the Type 1 quasar black hole mass function and Edding- 
ton ratio distribution, we merely note here that our ear- 
lier conclusions on the statistical properties of the broad 
line mass estimate errors are unchanged using a more 
flexible model for the Eddington ratio distribution. 

For the sake of brevity, in the remainder of this work 
we present the results obtained from the model that as- 
sumes that there is no luminosity-dependent bias, i.e., 
f3 = 0. We do this for easier comparison with other 
work, as almost all studies implicitly assume that there 
are no systematic trends in the mass estimate errors with 
luminosity. However, we also discuss what aspects of our 
results change when we allow the mean value of the error 
in the mass estimates to depend on luminosity. 

3.1. Evaluating the Fit Quality 

In order to evaluate whether our model provides a good 
fit to the dat a we use a technique called posterior pre- 
dictiv e check (|Rubinlll98il 119841 : iGelman. Meng. fc Sternl 
Il996h . For each random draw of the BHMF, the BHERF, 
and the values of a m i from their joint posterior probabil- 
ity distribution we generate a mock sample of virial mass 
estimates and luminosities, conditional on the SDSS 
quasar selection function. These mock samples are then 
compared to the actual data to check for consistency. 
This approach therefore includes our uncertainty in the 
BHMF and BHERF, as well as the randomness in gen- 
erating a single sample from the BHMF and BHERF. 

In Figures [2] and [3] we compare the histograms of the 
luminosities and virial mass estimates for our sample 
with the histograms of the mock samples generated from 
our MCMC output. Our estimated BHMF and BHERF 



QUASAR DEMOGRAPHICS II 



9 



10000 
1000 
100 
10 

1000 

S ioo 

5 io 
o 

S iooo 

E ioo 



1000 
100 
10 

I 




\ 



|z = 2.65 I 

5_ 



EEhD 



fl |z=1.60 I 



j |z = 3.20 | 




45 46 47 



45 46 47 48 



44 45 46 47 



45 46 47 48 

log XL X (2500A) [erg s" 1 ] 



Fig. 2. — Posterior predictive check comparing the actual dis- 
tribution of luminosities for our sample (solid histogram) with the 
set of distributions generated by our black hole mass and Edding- 
ton ratio functions (red squares with error bars). The error bars 
contain 68% of the posterior probability. The distributions of lu- 
minosities generated by our model are consistent with the observed 
distributions, showing that our model provides an acceptable fit. 
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Fig. 3. — Same as Figure[2] but for the virial mass estimates from 
IShen et al.l 1 1201 II ). The distributions of virial masses generated by 
our model are consistent with the observed distributions, showing 
that our model provides an acceptable fit. 



generate samples of the virial mass estimates and lumi- 
nosities that are consistent with our actual data, show- 
ing that our model provides an acceptable fit. The 
mock samples generated by the model that includes a 
luminosity-dependent bias term were also consistent with 
the data. 

3.2. Comoving Number Densities as a function of Mbh 

and z 

Our estimated BHMF and its uncertainty in various 
redshift bins is shown in Figure 2] and reported in Table 
[U The best-fit BHMF is defined as the posterior me- 
dian value of the number density as a function of Mbh, 
and the uncertainty in the BHMF is defined to be the 
region containing 68% of the posterior probability for 



TABLE 1 

Type 1 Quasar Black Hole Mass Functions 



z 


log M B h 


log <3>_ 


a log$ b 


log<J> + c 




(M e ) 


(Mpc~ 3 dex~ 




0.400 


8.000 


-5.00 


-4.64 


-4.21 


0.400 


8.100 


-5.03 


-4.67 


-4.30 


0.400 


8.200 


-5.06 


-4.72 


-4.39 


0.400 


8.300 


-5.09 


-4.78 


-4.50 


0.400 


8.400 


-5.14 


-4.87 


-4.62 


0.400 


8.500 


-5.22 


-4.97 


-4.76 


0.400 


8.600 


-5.32 


-5.12 


-4.92 


0.400 


8.700 


-5.47 


-5.29 


-5.09 


0.400 


8.800 


-5.65 


-5.49 


-5.31 


0.400 


8.900 


-5.87 


-5.72 


-5.54 


0.400 


9.000 


-6.12 


-5.97 


-5.80 


0.400 


9.100 


-6.40 


-6.25 


-6.08 


0.400 


9.200 


-6.70 


-6.54 


-6.37 


0.400 


9.300 


-7.02 


-6.84 


-6.67 


0.400 


9.400 
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Note. — The full tabic is available in the electronic 
version of the paper. Tabulated here are the results 
for the /3 — error model. Results for the (3 ^ error 
model arc similar to those presented in Paper I. 
a Lower boundary on the region containing 68% of the 
posterior probability for the BHMF, i.e., the 16 th per- 
centile of the posterior distribution. 
b Posterior median for the BHMF. 

c Upper boundary on the region containing 68% of 
the posterior probability for the BHMF, i.e., the 84 th 
percentile of the posterior distribution. 

the number density as a function of Mbh- Note that 
this represents the formal statistical uncertainty on the 
BHMF; in reality, the actual uncertainty is larger due to 
unaccounted for systematic errors. Also shown in Figure 
2] is the estimated 10% completeness limit and the flux- 
limited BHMF for quasars at i < 19.1 at z < 2.9 and 
i < 20.2 at z > 2.9, corresponding to the approximate 
flux limits of the SDSS DR7 quasar sample. 

It is apparent that our sample starts to become signifi- 
cantly incomplete at Mbh ^ 3x 10 8 M q , and any conclu- 
sions from the BHMF at these masses become reliant on 
our model for the BHMF and our assumptions for deriv- 
ing it from the mass estimates. In addition, extrapolation 
in our BHMF is unstable against small unidentified errors 
in both the SDSS selection function and the mass esti- 
mates beyond the region where we are reasonably com- 
plete. As discussed in K10, to zeroth-order the BHMF 
can be estimated as (I>m(Mbh) ~ ti(Mbh)/ s (Mbh), 
where u(Mbh) is the measured number density of Type 
1 quasars with mass Mbh, and s(Mbh) is the Type 
1 quasar completeness as a function of Mbh- The 
term s(Mbh) is calculated by averaging the SDSS se- 
lection function over the luminos i ty d istribution at fixed 
M BH (see Eq.[ll] in lKellv et all (|2009[ )). The estimated 
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Fig. 4. — The black hole mass function for broad line AGN derived from our sample. The blue region contains 68% of the probability for 
the BHMF, and the solid line running through it denotes the best-fit BHMF, defined to be the posterior median. The orange shaded region 
containes 68% of the probability on the BHMF for Type 1 quasars above the SDSS flux limits, and the dashed line denotes the posterior 
median. The vertical solid line denotes the mass below which the completeness in the BHMF drops to < 10% for the flux-limited SDSS 
sample, and thus mass bins below this marker are highly incomplete and the extrapolation becomes highly uncertain. We do not find any 
evidence for a turnover in the Type 1 quasar BHMF above the 10% completeness limit in Mbh- 



BHMF is highly sensitive to both the measured number 
of sources in a mass bin and the value of s(Mbh) when 
s(Mbh) is very small, and thus is unstable against even 
small errors in h(Mbh) and the selection function. 

The Type 1 quasars in the most incomplete bins tend 
to be those with the lowest S/N spectra. Because the 
FWH M measurement ca n be biased in low S/N spec- 
tra (|Dennev et al.l l2009h . the measured number densi- 
ties of the most incomplete mass bins may be incorrect 
due to biases in the mass estimates derived from FWHM 
measurements, which would scatter objects into incorrect 
mass bins. In addition, the SDSS quasar selection func - 
tion is derived via simulation by iRichards et al.l (|2006l ). 
Aspects of this simulation which are not representative of 
the quasar population (e.g., lack of a host galaxy compo- 
nent) can introduce errors into the quasar selection func- 
tion, which are manifested as errors in s(Mbh)- There- 
fore, in order to limit the impact of systematics, we focus 
our discussion on the regions of parameter space where 
the SDSS DR7 is > 10% complete. 

In general, the BHMF falls off approximately as a 
power- law toward higher mass Mbh, with the BHMF 
being steeper at higher mass. The comoving number 
density of SMBHs in Type 1 quasars continues to in- 
crease all the way down to the masses corresponding to 
the fa 10% completeness limits, and thus we do not ob- 
serve a peak in the BHMF from the SDSS. Instead, we 
constrain the peak in the Type 1 quasar BHMF to occur 



at M BH <2x 10 8 M q . 

In Figure [5] we show the evolution in the Type 1 quasar 
BHMF at four different values of Mbh- The number den- 
sities of SMBHs in Type 1 quasars with M B h < 10 9 M Q 
fall off more steeply toward higher redshift at z > 2, 
while the number densities of Mbh ^ 10 9 M Q fall off 
more steeply with decreasing redshift at z < 2. These 
trends imply that more massive black holes are more 
likely to be observed at higher z in Type 1 quasars than 
less massive ones. These trends have been called 'down- 
sizing', and have been obse rved in previous work (e.g., 
IVestergaard fc Osmerll2009L K10, Paper I). In addition, 
there is evidence for a discontinuity in the number den- 
sities across the the redshift where we switch from H/3 to 
Mgn when calculating the mass estimates. This discon- 
tinuity may be partly driven by systematic differences in 
the H/3- and Mgn-based mass estimates. However, it is 
also likely driven at least in part by the contribution from 
a host galaxy component to the nuclear emission which 
boosts the number of Type 1 quasars above the flux limit 
at z S, 0-8, but is not modeled in the selection function 
([Richards et alJl2~006| ). Considering this, it is likely that 
the number densities at z < 0.8 are overestimated, espe- 
cially at z = 0.4. 

In Figure [6] we compare the mass functions derived 
from both models for the error distribution in the virial 
mass estimates at three representative redshifts. There 
is little difference in the BHMFs derived from our mod- 
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Fig. 5. — Evolution in the comoving number densities of Type 
1 quasars at four different values of Mbh- The bars contain 68% 
of the posterior probability. The green bars denote bins with H/3- 
based virial mass estimates, the blue bars denote bins with Mgll- 
based mass estimates, and the red bars denote bins with ClV-based 
mass estimates. All of the bins show a peak in the number density 
around z ~ 2, with the higher mass bins falling off more steeply to- 
ward lower redshift and less steeply toward higher redshift, relative 
to the peak. The trend is commonly referred to as 'downsizing'. 



els with and without a luminosity dependent bias in the 
mass estimates, with the exception that the BHMF de- 
rived from Mgn and Civ is shifted toward lower values 
of Mbh by w 0.2 dex for the model which includes a 
luminosity dependent bias (i.e., j3 ^ 0). In addition 
the uncertainties in the BHMF arc larger when we al- 
low for a luminosity-dependent bias. We also note that 
the BHMF for the (3 ^ error model is consistent with 
that presented in Paper I, but with larger uncertainties 
on account of the more flexible Eddington ratio distribu- 
tion model. 

While the derived BHMF does not depend strongly on 
the existence of a luminosity-dependent bias, the flux- 
limited BHMF does depend strongly on the value of f3. 
This is because the segment of the SMBH population 
that is probed within a certain luminosity range depends 
strongly on how the errors in the mass estimates scale 
with luminosity. However, this is not a concern as the 
flux-limited BHMF is not of interest in this work. How- 
ever, it would be a concern, for instance, when using 
virial black hole mass estimates from flux-limited sam- 
ples to investigate models for accretion flows or evolution 
of the BH-bulge sc aling relation at high-redshift (e.g., 
iShen fc Kell^[20inh . 

3.3. How Massive can Black Holes Become?' 



Following K10, we compute the probability distribu- 
tion of Mgff(z), the most massive SMBH from a popu- 
lation of Type 1 quasars drawn from the BHMF within 
each redshift bin. In other words, this quantity may be 
thought of as the most massive SMBH that would have 
been observed within each redshift bin in an all-sky sur- 
vey of Type 1 quasars without a flux limit. The proba- 
bility distribution for the maximum value of Mbh gen- 
erated by a random draw from the BHMF is calculated 



log M'b^ 



p(mBH\Tr, fJ-, E) dmBH 



(19) 



where p(mBH\^, (J-, E) is obtained by marginalizing 
Equation ([2]) over luminosity. In order to incorporate our 
uncertainty on the derived BHMF, we calculate Equation 
(fTTjj) for each value of (N, tt, (i, E) returned by our MCMC 
sampler and average the results. 

Our derived constraints on M^J^z) are shown in Fig- 
ure [7] There are no obvious trends between Mgff(z) 
and redshift, although there may be a slight trend for 
Mgjf(z) to be larger at higher z. This is consistent with 
'downsizing', but may also be due to larger systematic 
errors in Civ-based mass estimates. The results for the 
luminosity-dependent bias model are very similar, but 
the uncertainties are ~ 30% larger. 



As mentioned in § l2.2l our derived constrained on Mgjy x 
are sensitive to the prior placed on Em,/c- This is because 
large value of E^fe imply tails in the BHMF that extend 
to larger value of Mbh, and therefore larger values of 
Mjjjf. If we constrain T*M,k < 2 instead of Em,Zc < 1, 
then the upper boundary of the error bars on M^ff ex- 
tend to M B h ~ 10 12 -10 13 M©. Ideally the prior should 
not have a strong influence on this quantity, but unfor- 
tunately the upper flux limit for the SDSS results in the 
extreme high mass tail of the BHMF being poorly con- 
strained; because of the upper flux limit, there is nothing 
in the data from prohibiting a small population of Type 
1 quasars hosting SMBHs with, say, Mbh ~ 10 12 M©. 
Considering this, we can formally only place a lower 
bound of Mgff ~ 1O 1O M0 on the maximum mass in the 
observable universe of a SMBH in a Type 1 . This is con- 
sistent with the range 2 x 10 10 M Q < Mf jf < 5 x 10 10 M Q 
at z > 1 that K10 estimated from their sample based on 
the SDSS DR3. 

3.4. Completeness in Black Hole Mass 

In Figure [5] we show the 10% completeness limit in 
Mbh as a function of z for a Type 1 quasar survey with 
a flux limit of i < 20 and i < 24. In calculating these 
completeness limits we assume a step function down to 
the flux limit, i.e., we do not assume the SDSS selection 
function. The z < 20 flux limit roughly corresponds to 
the SDSS flux limit, while the z < 24 flux limit is similar 
to the limiting magnit ude for the Pan-ST ARRS Medium 
Deep Survey Fields (|Saglia et al.1 12012?) as well as the 
spectroscopic samples fr om COSMOS (e.g., iLilly et all 
120071 iTrump et all [2007m . As with other quantities de- 
rived in this section, the results from the model with 
a luminosity-dependent bias in the mass estimates are 
similar but with larger uncertainties. Samples with a 
limiting magnitude of i — 20 start to become strongly 
incomplete at Mbh ~ 10 8 M Q by z ~ 1, increasing 
to Mbh ~ 7 x 10 8 M Q by z ~ 2. However, Type 1 
quasar samples with a limiting magnitude of i — 24 are 
able to go an order of magnitude 'deeper' in Mbh, be- 
coming incomplete at Mbh ~ 10 7 M Q by z ~ 1 and 
M BH ~ 7 x 1O 7 M by z ~ 2. 

The contribution from the host-galaxy is neglected in 
these calculations, but would likely become important 
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Fig. 6. — Comparison of the BHMF derived from the model with a luminosity dependent bias (fi 7^ 0, diagonal line-filled region) and 
without a luminosity dependent bias (/3 = 0, cyan region). For simplicity, we only show one redshift for each emission line used to calculate 
the virial mass estimates; the omitted bins exhibit similar differences. There is no apparent difference in the BHMFs derived from the H/3 
line, while the BHMFs derived from the Mgll and CIV lines under the model with /3 ^ are shifted lower in Mbh by 0.2 dex. This is a 
rather small difference and in general our scientific conclusions do not change when we allow a luminosity-dependent bias. 
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Fig. 7. — Constraints on the maximum mass that could be ob- 
served in a broad-line quasar as a function of redshift. Symbols 
are as in Figure [5] The maximum value of Mbh for quasars im- 
plied by our BHMF is Afgjf ~4x 10 10 M Q , although values of 
I0 10 < Mgff /Mq < I0 11 are consistent within the uncertainties. 
However, as discussed in the text the values of Mgr^ we derive are 
highly sensitive to our adopted prior constraints, and the values 
shown should be treated as lower limits. 



near i ~ 24, at least at lower redshift. The host-galaxy 
contribution would affect these calculations in two ways. 
First, the host galaxy leads to an increase in nuclear emis- 



FlG. 8. — The values of Mbh at which a Type 1 quasar survey 
becomes only 10% complete under a limiting magnitude of i < 20 
(red solid line) and i < 24 (blue dashed line), as a function of 
redshift. The data points denote the posterior median values, and 
the error bars contain 68% of the posterior probability. The SDSS 
becomes highly incomplete below Mbh ~ 3 x 1O 8 M0. 



sion, possibly moving the nuclear flux above the limiting 
magnitude, allowing one to detect intrinsically fainter 
Type 1 quasars. And second, the host galaxy dilutes the 
AGN emission, making it harder to identify whether the 
galaxy nucleus hosts a Type 1 quasar. The exact details 
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of these effects will depend on the distribution of Type 
1 quasar host galaxy luminosities and morphologies at 
fixed Mbh-, a s well as the Type 1 quasar identification 
algorithm of a particular survey. As such, it is unclear 
how the host galaxy would affect the Mbh completeness 
limits, and our estimated limits for i < 24 should merely 
be viewed as suggestive. 

4. THE EDDINGTON RATIO DISTRIBUTION 

4.1. Number Densities 

Our estimated Type 1 quasar black hole Eddington 
ratio function (BHERF) is shown in Figure [9] and re- 
ported in Table[2j Also shown is the flux-limited BHERF 
and the value of L/LiEdd below which our sample be- 
comes < 10% complete. As in Paper I, the BHERF was 
calculated from the model BHMF and luminosity func- 
tion assuming a bolometric correction to ^L ;y (2500A) of 
C2500 = 5. In general, the SDSS quasar sample is incom- 
plete at L/LEdd % 0.07. The comoving number densi- 
ties of SMBHs in Type 1 quasars increase toward lower 
Eddington ratio, a trend which continues beyond the in- 
completeness limit. The only exceptions are the z = 3.75 
and z — 4.25 bins, which display evidence for a peak in 
the BHERF at L/L Ed d ~ 0.3. However, the BHERFs 
in these bins are derived from the Civ line and may be 
subject to syste matics resulting from an outflowing win d 
component (e.g. JShen et al.ll2008t [Richards et al.ll2011h . 
so it is unclear if this peak is real. Indeed, the redshift 
bins before and after these two do not show any evidence 
for a peak in the BHERF, although they also employ the 
Civ line. In addition, we note that because we adopted a 
more flexible model for the Eddington ratio distribution, 
our constraint on the redshift evolution of the mean Ed- 
dington ratio is much poorer than that in Paper I, but is 
generally consistent with the trend found in Paper I. 

Figure [10] shows the evolution in the comoving number 
densities at four different values of L/LEdd- The number 
densities at L/L E dd ^ 0.05 are roughly constant at z < 
2, implying that there is no 'downsizing' in Eddington 
ratio at these redshifts. However, the number densities of 
Type 1 quasars radiating at L/LEdd ^5 0.1 decline rapidly 
at z > 3 while the number densities of Type 1 quasars 
radiating at L/LEdd <; 0.5 are similar at z < 2 and z > 3. 
This may be evidence for downsizing in Eddington ratio, 
in which the number densities of Type 1 quasars at low 
Eddington ratios increases rapidly from z ~ 4-5 to z ~ 3. 
However, we caution that this trend is primarily driven 
by the z — 3.75 and z — 4.25 redshift bins, which are 
the only ones that show a peak in the BHERF above the 
SDSS completeness limit. 

There is an apparent discontinuity in the number den- 
sities across z ~ 2, with the number densities at z ~ 2.15 
being ~ 1-2 orders of magnitude larger than that at 
z ~ 1.8. While the uncertainties on the number den- 
sities at 2 < z < 3 are large, making this only a l-2cx 
effect, there are a couple of other issues with this redshift 
range that are worth commenting on. For one, z ~ 2 
marks the transition between mass estimates calculated 
from Mgn and Civ, so there may be systematic differ- 
ences among these lines. However, we do not see the 
same effect in the evolution of the BHMF, and the num- 
ber densities of the virial mass estimates (i.e., the binned 
BHMF) typically show continuity between the redshift 
bins where the mass estimates switch emission lines (Pa- 



TABLE 2 

Type 1 Quasar Black Hole Eddington Ratio 
Functions 



z \ogL/L Edd log$_ a log<I>o b log*+ c 

(Mpc- 3 dex- 1 ) 

0.400 -1.50 -4.83 -4.45 -3.89 

0.400 -1.40 -4.84 -4.52 -4.00 

0.400 -1.30 -4.86 -4.60 -4.15 

0.400 -1.20 -4.91 -4.69 -4.31 

0.400 -1.10 -5.00 -4.83 -4.48 

0.400 -1.00 -5.13 -4.98 -4.68 

0.400 -0.90 -5.31 -5.18 -4.91 

0.400 -0.80 -5.54 -5.41 -5.17 

0.400 -0.70 -5.83 -5.70 -5.45 

0.400 -0.60 -6.17 -6.02 -5.76 

0.400 -0.50 -6.57 -6.40 -6.11 

0.400 -0.40 -7.03 -6.81 -6.49 

0.400 -0.30 -7.54 -7.25 -6.90 

0.400 -0.20 -8.09 -7.72 -7.31 

0.400 -0.10 -8.65 -8.18 -7.70 

0.400 0.000 -9.22 -8.63 -8.08 

0.400 0.100 -9.79 -9.08 -8.46 

0.400 0.200 -10.3 -9.49 -8.77 

0.400 0.300 -10.8 -9.88 -9.07 

0.400 0.400 -11.3 -10.2 -9.37 

0.400 0.500 -11.9 -10.6 -9.61 

0.400 0.600 -12.5 -10.9 -9.85 

0.400 0.700 -13.1 -11.3 -10.0 

0.400 0.800 -13.6 -11.6 -10.3 

0.400 0.900 -14.1 -12.0 -10.5 

0.400 1.000 -14.6 -12.3 -10.7 

Note. — The full tabic is available in the electronic 
version of the paper. Tabulated here arc the results for 
the /3 — error model. Results for the /3 ^ error model 
arc similar to those presented in Paper I. 

a Lower boundary on the region containing 68% of the 
posterior probability for the BHERF, i.e., the 16* % per- 
centile of the posterior distribution. 
b Posterior median for the BHERF. 

c Upper boundary on the region containing 68% of the 
posterior probability for the BHERF, i.e., the 84 th % per- 
centile of the posterior distribution. 

per I) , suggesting that the use of different emission lines 
is not the dominant reason for this discontinuity. In- 
stead, the apparent discontinuity in number density may 
be primarily due to systematic errors in our incomplete- 
ness correction. The three redshift bins z = 2.15,2.65, 
and 3.20 correspond to redshifts where quasars colors 
are similar to star colors, making the SDSS quasar color 
selection incomplete at these redshifts. The complete- 
ness can be as low as ~ 5% in this redshift range. The 
color distribution of simulated quasars does not perfectly 
match the observed color distribution at these redshift 
((Richards et al.|[2006l) , suggesting that there may be sys- 
tematic uncertainties in the estimated selection function. 
In addition, the completeness of the selection algorithm 
at these redshifts depends on the quasar optical/UV 
SED, which in turn has bee n found to depend on L/Lsdd 
(e.g.. iBonning et all l2007t ). However, any dependence 
of this completeness on L/L Edd is not accounted for in 
the SDSS selection function, introducing systematic er- 
ror if the quasar colors do depend on L/LEdd- Similarly, 
the Type 1 quasar bolomet ric correction also depends 
on L/LEdd and Mbh (e.g.. IVasudevan fc Fabianl 120071 : 
Kelly et al.ll2008t IVasudevan fc Fabianl 12009) . Our use 
of a constant bolometric correction may exasperate the 
systematics over this redshift range. Because the esti- 
mated number densities are unstable in bins of L/LEdd 
that are significantly incomplete (see £|3.2[) . there is likely 
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Fig. 9. — The black hole Eddington ratio function for broad line AGN derived from our sample. The labeling is the same as in Figure 
[9] There is a broad range in L/L^dd f° r Type 1 quasars, and there is no evidence for a turn-over in the BHERF down to L/L^^ ~ 0.07, 
except for possibly in the redshift range 3.20 < z < 4.75. 
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Fig. 10. — Evolution in the comoving number densities of Type 
1 quasars at four different values of L/L^dd- The labeling is the 
same as in Figure [5] The number densities of Type 1 quasars are 
fairly constant in redshift below z < 2 for all Eddington ratio bins. 
However, the number densities of Type 1 quasars radiating at lower 
values of L/L^dd drop off more steeply toward higher redshift at 
z > 2. 



a significant additional systematic uncertainty that is not 
reflected in the error bars on the number densities at 
2 < z < 3.2. Therefore, we do not consider the apparent 
discontinuity in the number densities at z ~ 2.15 to be 
real. 

Unlike with the BHMF, the BHERF is noticeably dif- 



ferent when we allow for a luminosity-dependent bias. 
In Figure [11] we compare the BHERF derived from the 
models with and without a luminosity-dependent bias 
for three representative redshift bins. This difference is 
negligible for the H/3-based mass estimates, for which 
/3 « 0. However, in the redshift bins where Mgn and 
Civ mass estimates are used the model that includes the 
luminosity-dependent bias leads to BHERFs which are 
more uncertain and fall off flatter toward higher L/ LEdd, 
implying a larger number of Type 1 quasars radiating 
near the Eddington limit. In addition, while the number 
densities at the high L/Lsdd end of the BHERF depend 
on the value of /3, the evolution in the number densities 
at fixed L I Lsdd is not as affected by the value of (3. The 
evolution results obtained for the luminosity-dependent 
bias model are not significantly different from those ob- 
tained for the j3 = model. 

4.2. Distribution of Eddington Ratio as a Function of 
Black Hole Mass 

In Figure [T2l we show the median and 95 th percentile 
of the distribution of Type 1 quasar Eddington ratio at 
fixed Mbh', note that the 95 th percentile is the value of 
L/LEdd\MBH at which 95% of quasars would have be 
below this value. At most redshifts the median value 
of L/LEdd\MBH is below our incompleteness limits, so 
we limit our discussion to trends involving the 95 th per- 
centile of the L/ LEdd\MsH distribution; the median val- 
ues are simply shown for reference and in general exhibit 
similar trends as the 95 th percentile. 

At both high (z > 3.2) and low (z < 0.6) redshifts the 
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Fig. 11. — Same as Figure [6] but for the Eddington ratio function. Unlike for the BHMF, there is a more significant difference in the 
BHERFs derived under the two models for the virial mass estimate errors. In all cases the uncertainties on the BHERF are larger for the 
mode that includes a luminosity-dependent bias. In addition, the BHERF under the model with ^ is shifted toward higher values 
L/L^dd- This shift is small in the redshift bins using H/3, but increases to 0.5 dex for the redshift bins using Mgll and CIV. 



value of the 95 th percentile in the L/LEdd distribution 
is relatively independent of M BH . However, at redshifts 
0.8 < z < 2.65 the value of the 95 th percentile of the 
distribution of L/LEdd at fixed Mbh increases toward 
larger values of Mbh- This therefore implies that at 
0-8 J$ z < 2.65 Type 1 quasars with more massive black 
holes are more likely to be radiating near the Edding- 
ton limit. However, we note that this redshift range is 
dominated by mass estimates derived from Mgn, and it 
is possible that systematic effects with this line may be 
driving some of these results. That being said, the results 
in the last two redshift bins in this range are derived from 
the Civ line mass estimates, which show the same trend 
of a decreasing fraction of Type 1 quasars radiating near 
the Eddington limit with decreasing Mbh- The results 
from this section are similar for the model that includes 
a luminosity-dependent bias, with the exceptions that 
the uncertainties in the 95 th percentiles were larger, and 
that the value of the 95 th percentile does not increase as 
steeply with Mbh over 0.8 < z < 2.65. 

These results are in contrast to the conclusions reached 
bv lSteinhardt fc El vis (2010a). These authors used virial 
mass estimates fr om the SDSS DR5 quasar sample of 
IShen et all (|2008[ ) to argue that the most massive Type 
1 quasars fall short of the Eddington limit, forming 
what they called a 's ub-Eddington boundary' (b ut see 
IRafiee fc Hall 20llah . iSteinhardt fc Elvisl l|2010ah quan- 
tified this trend through the 95 th percentile of lumi- 
nosities above the peak in the distribution of estimated 
bolometeric luminosities for the sample, finding that the 



95 percentile in luminosity increased slower with Mbh 
than would be expected from a constant value of L/LEdd- 
We consider the likely reason for their different con- 
clusion to be uncorrected incompleteness. While it is 
true that for a flux-limited sample the 95 th percentile in 
L/LEdd at fixed Mbh mass does increase with decreasing 
Mbh, this is not necessarily true of the Type 1 quasar 
population as a whole. The reason for this is because 
the 95 th percentile in L/LEdd at fixed Mbh for a flux- 
limited sample is an overestimate of the true 95 th per- 
centile due to the loss of the faint end of the population. 
For smaller values of Mbh one loses a larger fraction 
of the lower L/LEdd part of the population, increasing 
the bias in the 95 th percentile inferred from the distribu- 
tion of Type 1 quasars that are actually bright enough 
to be detected. Indeed, for the least massive black holes 
in a flux-limited sample one cannot even detect the 95 th 
percentile of the Eddington ratio distribution, as it falls 
below the flux limit. This leads to a spurious increase in 
the 95 th percentile of the L/LEdd distribution with de- 
creasing Mbh- In addition, the errors of virial BH mass 
estimates stretch the distribution in the mass-luminosity 
plane along the mass direction, causing artificial flatten- 
ing and an apparent "tilt" away from the Eddington limit 
towards higher virial BH masses (see Fig. 8 of Paper I), 
which could be incorrectly recognized as visual evidence 
for a sub-Eddington boundary. However, we note that 
this effect is not as strong as the bias caused by incom- 
pleteness. 

In this work we have corrected for incompleteness in 



16 



KELLY & SHEN 




0.001 1 1 1 

8 8.5 9 9.5 10 10.5 8.5 9 9.5 10 10.5 11 

log M BH /M @ 

Fig. 12. — Dependence of the location of the median (orange regions) and 95 th percentile (blue regions) in the Eddington ratio distribution 
at fixed Mgjj and redshift. The solid and dashed lines denotes the posterior median value, and the shaded region contains 68% of the 
posterior probability. The high Eddington ratio tail of the distribution does not change with Mbh at z < 0.6 and z > 3.2, but moves to 
higher values of L/L^dd with increasing Afgjy at 0.8 < z < 3.2. This unusual behavior may be due to unidentified systematics in the virial 
mass estimates. 



flux-limited samples and errors in virial BH mass esti- 
mates, assuming our statistical model and the SDSS se- 
lection function, and thus attempt to recover the true in- 
trinsic trends in the high L/LEdd tail with Mbh- Upon 
doing so, we do not find any evidence that higher mass 
Type 1 quasars are less likely to be observed at high 
L/LEdd, but in fact at certain redshifts may be more 
likely to have high L/LEdd- Therefore, the observed 
dirth of Type 1 quasars having both high M bh and high 
L/Lsdd is not caused by a change in the shape of the tail 
of the Eddington ratio distribution, but is instead caused 
by the rapidly decreasing number densities of Type 1 
quasars toward higher masses. Because the most mas- 
sive SMBHs in Type 1 quasars are rare by definition, 
we would not expect those that we do observe to also oc- 
cupy the rare high-L / L Edd region of the mass-luminosity 
plane. 

For full disclosure, at low redshift past the peak of 
quasar activity, the most massive SMBHs are probably 
accreting mostly at very low Eddington ratios, but are 
no longer shining in the form of Type 1 quasars. We 
do not find any statistical evidence of a sub-Eddington 
boundary in the mass-luminosity plane of type 1 quasars, 
and earlier claims of such a boundary are likely the result 
of uncorrected incompleteness. 

4.3. Do quasars obey the Eddington Limit? 

Similar to the calculation of the maximum mass per- 
formed in § 13.31 we can also calculate the probability 



distribution of the highest value of L/LEdd from a pop- 
ulation of Type 1 quasars drawn from our estimated 
BHERF. The maximum value of L/LEdd as a function of 
redshift is shown in Figure [T3] Our BHERF implies that 
the maximum value of L/LEdd for a Type 1 quasar is 
L/LEdd ~ 3, although values between L/LEdd — 1 and 
L/LEdd = 10 are well within the uncertainties. The re- 
sults from using the model with a luminosity-dependent 
bias suggests values of the maximum Eddington ratio a 
factor of ~ 2 higher, which is smaller than the uncertain- 
ties within individual redshift bins. The maximum values 
of L/LEdd are somewhat divergent at 2 < z < 3.2, having 
values near L/LEdd ~ 10, but, as discussed above, this 
range is subject to large statistical and systematic un- 
certainties. In addition, our results suggest that quasars 
at z ~ 4 may obtain a maximum value of L/LEdd that is 
slightly larger than quasars at z < 2, but the error bars 
are large. Considering this, our results are consistent 
with a maximum value of L/L Edd ~ 3 across all redshift 
bins. In addition, there is additional systematic uncer- 
tainty in the maximum value of L/LEdd caused by our 
assumption of a constant bolometeric correction. Taking 
this into account, we do not find any evidence that Type 
1 quasars violate the Eddington limit by more than a 
factor of a few. 

4.4. Completeness in Eddington Ratio 

Similar to the calculation performed in § 13.41 in Fig- 
ure [14] we show the Eddington ratio at which a survey 
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Fig. 13. — Constraints on the maximum Eddington ratio that 
could be observed in a broad-line quasar as a function of redshift. 
Symbols are as in Figure [5] There is no evidence for significantly 
super Eddington radiation, and the maximum value of L/LEdd f° r 
quasars implied by our BHERF is L/LEdd ~ 3. 



Fig. 14. — The values of L/LEdd at which a Type 1 quasar 
survey becomes only 10% complete under a limiting magnitude of 
i < 20 and i < 24. The labeling is the same as in Figure [8] The 
SDSS becomes highly incomplete below L/LEdd ~ 0.07. 



becomes only 10% complete as a function of redshift for 
two limiting z-band magnitudes. The results are similar 
for the model with a luminosity-dependent bias. Surveys 
with a limiting magnitude of i = 20 are largely incom- 
plete at L/LEdd i$ 0.07, while surveys with a limiting 
magnitude of i = 24 start to become severely incomplete 
at L/LEdd ^$ 0.01. The increase in the completeness 
limits at 2 < z < 3.2 mirror other anomalous trends in 
L/LEdd, and the uncertainties are large. As discussed 
earlier in this work, the SDSS color selection algorithm 
has difficulty distinguishing quasars from stars in this 
redshift range, and this may introduce systematic trends 
with L/LEdd- Therefore, we do not consider the 'bump' 
in the completeness limit at 2 < z < 3.2 to be real. These 
derived completeness limits are broadly consistent with 
the lack of Type 1 qu asars having L/Lem < 0.01 ob- 
served in COSMOS bv lTrump etaLl (|2011h . whose sam- 
ple starts to become incomplete at i > 23. 

5. DISCUSSION 

5.1. Black Hole Growth Time Scales: Evidence for 
Self-regulated Growth of SMBHs? 

In Figure [15] we show the typical growth time as a 
function of redshift for SMBHs in Type 1 quasars with 
Mbh/M & = 5xl0 8 ,10 9 ,5xl0 9 , and 10 10 . These growth 
times are cal culated assuming a radiative efficiency o f 
e r = 0.1 (e.g. JYu fc Tremaindl200llDavis fc Laorll2011h . 
and a seed mass of Mbh = 10 6 Af Q at z = 10. The 
implicit assumption in these calculations is that the 
SMBHs in each mass bin accrete at a time-averaged 
rate relative to Eddington that is equal to the popu- 
lation average of Type 1 quasars in that bin. While 
this is not true for every Type 1 quasar in each bin, 
it should be a good approximation for a representa- 
tive quasar in each bin; hence the association of these 
time scales with the 'typical' growth time for a SMBH 
in that bin. For comparison, we also show the age of 
the universe as a function of z. The results were simi- 
lar for the model with a luminosity-dependent bias, but 
the error bars on the growth times were larger. The 
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Fig. 15. — Typical time needed to grow a black hole from a 
seed mass of M|j c ^ d = 1O 6 M to four different final masses as a 
function of redshift, assuming a radiative efficiency of e r = 0.1 
and the typical Eddington ratio for that mass bin. The solid line 
marks the age of the universe as a function of z, the dashed line 
marks the growth time assuming L/LEdd = b an d the rest of the 
labeling is the same as that in FigurcO Although the uncertainties 
are large, the growth times for SMBHs in Type 1 quasars having 
Meh J 5 X 10 s Mq are comparable to or greater than the age of 
the universe, suggesting an earlier phase of accelerated obscured 
growth. The growth times at z < 0.6 are also longer than the age 
of the universe, possibly reflecting a shift to fueling by mass loss 
from evolved stellar populations. 



value of e r — 0.1 corresponds to the radiative efficiency 
of a moderately spinning black hole, and the value of 
Mbh = lO 6 Af for the seed black holes corresponds to 
the largest values predicted by models for SMBH seed 
formation (for a review, see IVolonteril[2rjl0l ). In partic- 
ular, SMBH seed masses of Mbh = 10 6 M Q represent 
the high-end of the distribution predicted from models 
wher e SMBH seeds form from the direct collapse of gas 
(e.g.. IVolonteri et al.ll2008t IVolonteri fc Begelma3 l2010t 
iNataraian fc Volonterill2012ft 
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For many of the mass and redshift bins at z > 2 
the growth times are comparable to or longer than 
the age of the universe, although the uncertainties are 
large for several of the bins. If the remnants of Pop 
III stars are the SMBH seeds then values of the seed 
masses are M R „ ~ 1 00M Q (e.g.. lFrver fc Kalogerall2001l : 
iMadau fc~ Rccs 200l]) and the growth times will be a fac- 
tor of ~ 2 longer. In addition, if the black holes are 
rapidly spinning, as expected from continuous accretion 
in a disk geometry, then their growth time will be a 
factor of ~ 3 longer due to the increased radiative ef- 
ficiency. However, we have neglected mergers of black 
holes in these growth times, which may make a non- 
negli gible contribution t o shortening the growth times 
(e.g-. ISesana et alj|2007t iTanaka fe Haimanl 12001 . Sim- 
ilarly, if most quasars at high z host SMBHs that are 
not spinning the radiative efficiency is e r < 0.1 and the 
growth times will shorter than what we estimate. 

Because our estimated typical SMBH growth time 
scales are comparable to, or longer than, the age of the 
universe, they imply that many SMBHs underwent an 
earlier phase of accelerated growth (e.g., with higher Ed- 
dington ratio) when we would not have observed them as 
Ty pe 1 quasars. Similar conclusions have been reached 
by iNetzer et al.l (|2007al) . K10, and iTrakhtenbrot et al.l 
(|2011l ). It is possible that we do not observe these early 
rapidly growing SMBHs as Type 1 quasars because they 
are obscured. Alternatively, it is also possible that the 
typical Eddington ratio is near unity for Type 1 quasars 
with masses representative of an earlier phase in their 
growth (e.g., Mbh J; 5 x 10 8 M q ). In this case, we would 
have missed these rapidly accreting Type 1 quasars be- 
cause they would have had smaller masses, and therefore 
would have fallen below our completeness limit. Both 
possibilities exist but cannot be distinguished by our 
data. 

Our results regarding the growth times of typical 
SMBHs in Type 1 quasars is in agreement with ex- 
pectations fro m models for self-r e gulated black hole 
grow t h (e.g., . Sanders et al.l 119881 : ISanders fc Mirabel 
119961 IHopkins etall l2006bf l. Within the framework of 
these models, the black hole is enshrouded in material 
and undergoes Eddington limited obscured growth. The 
broad line quasar phase begins after the obscuring ma- 
terial disappears. This can happen because the black 
hole grows to the point that feedback energy from it 
becomes powerful enough to 'blow' the obscu ring ma- 
terial away (e.g., IHopkins et all 120051 l2006bl ). Alter- 
natively, the obscuring material may disappear because 
it becomes consumed by star-formation in the bulge. 
The Type 1 quasar phase is expected to persist until 
the accretion rate drops low enough to switch to a ra- 
diatively inefficient accretion flow (e.g.. IChurazoy et al.l 
l2tM ICacI 120051: iShen et al.ll2"007a! IHopkins et al.ll2009h . 
We note that during the blow-out phase there may still 
be geometry-dependent obscuration caused by a dusty 
torus, a s invoked by the traditional AGN unif ication 
models (|Antonuccil fl99llUrrv fc Padovanilll995l ). Our 
results are therefore consistent with a variety of models 
for growing the most massive SMBHs at z > 2, so long 
as SMBHs exhibited enhanced accretion rates relative to 
Eddington below M B h ~5x 10 8 M©. 



5.2. Comparison with other Work: Implications for the 
Fueling and Growth of Supermassive Black Holes 

The results that we have obtained in this work with 
regard to AGN demographics provide an important win- 
dow into SMBH growth and AGN fueling, especially 
when interpreted within the greater context of ear- 
lier observational and theoretical work on AGN num- 
ber densities, clustering, and host galaxy properties. 
In particular, a picture is emerging where the fuel- 
ing mechanisms of SMBH growth are complex and 
varied. In this section we discuss an interpretation 
of our results that incorporates other recent observa- 
tional results; many of the ideas discussed here are 
treat ed in greater detail in various theoretical papers 
(e.g.. IHopkins fc Hernaulstl I2006I: IHopkins et all [2 006b; 

Bower et al.l l2006t iMonaco et al.l 



Croton et ; all 120061; 
20071; IHopkins et al 



Fanidakis et al.ll201 



20081; IHopkins fc Hernquistl 120091; 



Due to incompleteness, our sample probes the Type 
1 quasar population with Mbh ^ 3 x 10 8 M Q . Lo- 
cally, SMBHs in this mass range are predominately 
found in early- type (i.e., elliptical and SO galaxi es) galax- 
ies a s inferr ed from the local SMBH BHMF (|Yu fc Lul 
12004 I2008D . This therefore suggests that our study 
is dominated by SMBHs that currently live in ellipti- 
cal and SO host galaxies. Indeed, many studies find 
that the host galaxies of optically-id entified AGN are 
dominated by early type galaxies (e.g.,[Ka uffman n et all 
l2003UZakamska et alJl2006l : lKotilainen et al.ll2007| ). Cur- 
rently, the favored mechanism for cre ating elliptical 
galaxies is through galaxy mergers (e.g.. ISpringel et all 
20051: iBoumaud et all 120051 : iBovlan-Kolchin et al.l 120061 : 
Cox et al.ll2006l) . suggesting that the Type 1 quasars in 



our sample have at least in part been fueled by merg- 
ers. This is supported by observations which find evi- 
dence of past i nteractions in many local elliptical host s 
of quasars (e.g.. lBahcall et allll997l:lBennert et~aTll2008ft . 
However, we also note that some recent theoretical work 
suggests that spheroid-dominated galaxies m ay also form 
from instabilities in gas-rich disks at h igh-z (|Dekel et all 
2009); furthermore. iBournaud et"all (|2011l) has argued 
that AGN at high- z may form from similar processes. 
This said, the fraction of galaxies at high-z that contain 
the massive clumps necessary to form sphero ids and AGN 
via s ecular processes is potentially low (jWuvts et all 
l2012f ). and at best unclear. 

As discussed above, the typical growth times that we 
find for M BH >5x 1O 8 M SMBHs in Type 1 quasars 
at z > 2 are of order the age of the universe or longer, 
implying that these SMBHs experienced an earlier phase 
of accelerated growth. This is expected from models for 
self-regulated SMBH growth. Such models posit a mas- 
sive fueling event, such as a major merger, which directs 
larger amounts of gas toward the nuclear region, fueling 
obscured Eddington-limited SMBH growth. Eventually, 
the SMBH's growth is quenched either due to AGN feed- 
back or the consumption of gas due to star formation, 
revealing the SMBH as a Type 1 quasar and linking the 
SMBH's final mass with properties of the stellar bulge. 
Indeed, it is likely that most of our Type 1 quasars cur- 
rently reside in early type g alaxies, which show t he tight- 
est M B H-o-* relation (e.g.. iGiiltekin et alll2009l ). More- 
over, during the Type 1 quasar phase the AGN lightcurve 
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exhibits a power-law like decay, either due to the decrease 
in the fueli ng rate due to evolution of a feedback-driven 
blast wave (|Hopkins fe Hernquistll2006l ) or due to viscous 
evolution of an a ccretion d i sk resulting from a quenche d 
fuel supply (e.g., lYu et all [20051 : iKmg fc Pringldl2007h . 
If the AGN lightcurve is decaying during the Type 1 
quasar phase, then we would expect a broad range of 
Eddington ratios with the number densities increasing 
monotonically toward lower values of L/LEdd- This is 
consistent with our estimated BHERF, although we note 
that predictions for this type of BHERF are not unique 
to self-regulated black hole growth models. 

At redshifts z < 0.8 our calculated black hole growth 
times are a factor of ~ 2-3 longer than the age of the 
universe for SMBHs with M BH > 5 x 1O 8 M , reflecting 
the smaller values of L/LEd d at these redshifts. Simi- 
lar results were obtaine d by iHeckman et all ([20041 ) and 
iKauffmann fe Heckmanl (|2009l ) using a sample of low- z 
Type 2 AGN. These low-z AGN likely experienced an 
earlier stage of rapid growth, and instead we are currently 
witnessing a reignition of weaker AGN activity reflected 
in their lower accretion rates and longer growth times. 
Because the Type 1 quasars probed by our study most 
likely live in elliptical galaxies, this low-order AGN ac- 
tivity is probably fueled by stellar mass loss from evolved 
stars. Although there is little overlap in Mbh between 
the massive systems in our sample and the le s s mas sive 
ones in the sample of IKauffmann fe Heckmanl (|2009f ) . at 
Mbh ~3x 10 8 M q they derive an Eddington ratio dis- 
tribution which decreases monotonically toward higher 
L/LEdd, similar to us. 

In this work we have also found evidence for down- 
sizing in both Mbh and L/LEdd- The downsizing in 
Mbh is consiste nt with results from SMBHs in Type 
1 qua sars (e.g., IVestergaard et all l2008t lLabita et al.l 
I2009ai K1 Paper I) and for the entire SMB H popu- 
lation (e.g., iMarconi et al.ll2004t iMerlonil l200l ) . and re- 
flects the fact that the most massive SMBHs tended to 
experience active phases and the bulk of their growth be- 
fore less massive SMBHs. This, combined with the fact 
that AGN are observed to reside in dark matter halos 
of M h ~3x 10 12 M W at all redshifts (e. g . IShen et al 



20071J |Sh en et al.| 120091; [Ross et al.| [200§ |Hickox et all Uroum et al. l_'()lt>: Uahm et al. l_'UU!): Scljawu^ki et al. 
2009Q led |Hickox et al] ([20091) to argue that AGN ac- [20Tlt iKocevski et al.|[2012T ). with the later studies having 



internal processes and experiencing weaker AGN activ- 
ity. However, because our sample becomes significantly 
incomplete below L/LEdd ^ 0.07 it is unclear how sig- 
nificant this trend is. 

Based on our derived BHERF, we conclude that those 
Type 1 quasars that do radiate near the Eddington limit 
are extremely rare, suggesting that if Type 1 quasars 
do violate the Eddington limit they do so only for a very 
brief period of time. In such super-critical accretion flows 
the luminosity depends logarithmically on the accretion 
rate, so in principle these Type 1 quasars that do radiate 
at L/LEdd > 1 could be accreting at a significantly higher 
rate relative to Eddington, i.e., M/MEdd 3> 1. 

The most massive SMBHs in the observable Universe 
implied by our BHMF have masses 1O 1O M Q -1O 11 M0, 
consistent with earlier results obtained by K10. However, 
as this value is strongly dependent on our prior for the ex- 
treme high mass tail of the BHMF, a more realistic con- 
strain t is M™jf > 1O 1O M . Recently IMcConnell et all 
(|201lD detected two SMBHs in the local universe with 
Mbh ~ 1O 1O M0 determined from stellar dynamical 
modeling. Considering the different volumes probed by 
the SDSS DR7 and their search, their results are con- 
sistent with the limits on the maximum mass implied by 
our BHMF. In addition, the lower limits on the maximum 
masses that we infer are consistent with values predicted 
from mod els that assume the SMBH 's growth is self- 
regulated (Na taraian fc Trei stcr 2009) . Simulations that 
follow the growth of SMBHs in bright z ~ 6 quasars have 
also been able to grow SMB Hs to M B h ~ 2 x 1O 1O M 
by z ~ 2 ([Siiacki et al.ll2009D . 

Due to incompleteness our sample only probes SMBHs 
with Mbh ^ 3 x 10 8 M Q , and it is unlikely that our 
results and discussion extend to lower mass SMBHs. Es- 
timates of the local BHMF for all SMBHs imply that 
late types should dominate SMBH hos t galaxies below 
5 x 10 7 M Q (e.g., lYu fc Lul l200i 12001 . Furthermore, 
among galaxies that have been used to define the Mbh- 
a* rel ationship, late types ar e found at Mbh ^ 10 8 M© 
(e.g., fMcConn ell et al.ll2011[) . Many recent studies of 
X-ray selected AGN out to z ~ 2 find that they re- 
side in both disk- and bulge- dominant galaxies (e.g., 
Grogin et al. 2005; Gabor et al. 2009; Schawinski ct al. 
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tivity is triggered when a SMBH's host halo reaches 
Mh ~3x 10 12 M Q . Within this interpretation, the ob- 
served downsizing in Mbh is a reflection of the fact that 
the most massive SMBHs are those that reside in the 
most massive halos at high redshift, and these are the 
halos that reach a critical mass of Mh ~ 3 x 10 12 M, 
first. 

The downsizing in L/LEdd has a somewhat different 
form than that of Mbh ■ Unlike the number densities of 
Mbh, the number densities in bins of constant L/LEdd 
do not show any evolution at z < 2. Instead, we observed 
'downsizing' in L/LEdd in the sense that Type 1 quasars 
radiating at L/LEdd 0.1 arc significantly more rare at 
z ~ 4 compared to z < 2, while the number densities of 
Type 1 quasars at L/LEdd > 0.1 are similar for z ~ 2 
and z ~ 4. This may reflect a stronger contribution 
of the BHERF from Type 1 quasars at z ~ 2 that are 
further along the decaying part of the AGN lightcurve, 
and thus at a lower L/LEdd, as well as a stronger contri- 
bution of Type 1 quasars that are fueled through more 



a fainter flux limit and fin ding a higher f r action of AGN 
in disk galaxies. Indeed, IKocevski et al.l (|2012f ) found a 
correlation between X-ray luminosity and the fraction of 
AGN in spheroidal galaxies, suggesting that AGN with 
more massive SMBHs tend to be found in early types 
even at z ~ 2. The larger fraction of AGN with less mas- 
sive SMBHs in disk galaxies suggests that these SMBHs 
are grown through secular processes up to the point that 
we observe them, and it is unclear if they also experi- 
ence an earlier phase of accelerated growth. Moreover 
the weaker Mbh—c * relationship for late type galaxies 
(e.g.. i Grahaml 120081; iGiiltekin etaTI 120091; iGreene et all 
I2010at iKormendy et al.ll2011l ) suggests that the growth 
of these SMBHs is more weakly coupled to the evolution 
of the host galaxy, if any relationship exists at all. In 
contrast, the more massive SMBHs (M B h > 3 x 1O 8 M ) 
probed by our Type 1 quasar sample are grown through 
a process that also generates a spheroid and places the 
SMBH on the Mbh-o'* relation, with a major merger of 
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two gas rich galaxies being a likely candidate. 

Direct comparison of our results with models of SMBH 
fueling and growth is difficult as most models do not 
present a quantitative prediction of the BHMF and 
BHERF specifically for Type 1 quasars. Instead, many 
modelers predict the BHMF and BHERF for all SMBHs, 
or for active SMBHs. Generically, many models pre- 
dict a BHMF which increases toward lower mass down 
to Mbh S 3 x 10 8 Af©, which is in agreement with 
our estimated BHMF. In addition, many models predict 
downsizing in SMBH growth and AGN activity, in agree- 
ment with our empirical results. INataraian fc Volonteril 
(|2012t l compared BHMFs for Type 1 quasars with mod- 
els for merger driven SMBH growth that assumed that 
either SMBH seeds are the remnants of Pop-Ill stars 
or we re formed through direct collapse of pre-galactic 
disks (|Lodato fc Nataraianll200"7t ). They compared their 
model unobscured quasar BHMF to those estimated by 
K10, and in general they found that the Pop-Ill seeding 
model underpredicted the BHMF at M BH > 10 9 M Q . 
Their direct collapse seeding model provided a better 
fit at z < 4 although it overpredicts the number den- 
sities of the most massive SMBHs in Type 1 quasars at 
z < 2; however, the authors argue that the overpredic- 
tion is not a problem as it results from the fact that 
they did not implement depletion of the available gas to 
grow the SMBH at lower z. None of their models were 
able to match the K10 BHMF at z = 4.25. Because the 
K10 BHMF was derived from the SDSS DR3, which is 
a subset of our sample, it is not greatly different from 
the BHMF that we derive here. As such, th e conclu- 
sions reached by INataraian fc Volonteril (|2012f ) would be 
unchanged using our newer BHMF, although we note 
that comparison with our new BHMF implies that all of 
their models underpredict our derived number densities 
of SMBHs in Type 1 quasars with M BH ~ 3 x 10 8 M© at 
z = 1.25. 

5.3. Sources of Systematic Error 

While we have obtained a number of interesting re- 
sults, there are a few caveats regarding our approach 
that must be kept in mind, and we conclude this section 
with a discussion of them. There are three significant 
potential sources of systematic error in our approach, 
which we have touched on earlier: incorrect specification 
of the virial mass error distribution and bias, errors in 
the SDSS selection function, and errors in the bolome- 
teric correction. In this work we have used a simple 
constant bolometeric correction to the 2500A luminos- 
ity, and, strictly speaking, our derived BHERF should 
be viewed as the number density of Type 1 quasars radi- 
ating at a given ratio of optical luminosity to Eddington, 
scaled upward by a constant factor of five. Because the 
bolometeric correction likely depends on both Mbh and 
L/LEdd this will create systematic errors. It is unclear 
how these systematics in the bolometeric correction af- 
fect our results on the BHERF. In principle it is possible 
to incorporate a variable bolometeric correction into our 
Bayesian model, but unfortunately there is considerable 
uncertainty in the distribution of bolometeric corrections 
as a function of Mbh and L/LEdd- 

As discussed in § 13.21 and § 14.11 systematic errors on 
the selection function can have a significant effect on our 
results in highly incomplete regions of parameter space. 



In general we have limited our analysis to masses or Ed- 
dington ratios in which we are not highly incomplete (i.e., 
completeness > 10%), so we do not expect small errors 
in the selection function to have a significant effect on 
our results. However, this is not true in the redshift 
bins corresponding to z = (2.15, 2.65, 3.20). In these bins 
the SDSS selection algorithm has difficulty distinguishing 
quasars from stars, and as a result these bins are highly 
incomplete. Moreover, the distribution of quasar colors 
obtained from the simulations used to estimate the selec- 
tion function do not perfectly match the observed distri- 
butions, and the selection algorithm is color-dependent 
at these redshifts. We do not include a color-dependence 
in the selection function, and, because the quasar SED 
depends on Mbh and L/LEdd-, this likely introduces sys- 
tematic errors into our incompleteness correction. Be- 
cause of this the BHMF and BHERF derived at these 
redshifts should be interpreted with caution. 

The most important possible source of systematic error 
in our approach arises from incorrectly specifying the er- 
ror distribution of the virial mass estimates. This is par- 
ticularly a concern for Mgn and Civ, which are not as 
well studied with respect to the reverberation mapping 
database as H/3. In addition, Civ emission is thought 
to at least partly arise from an accretion disk wind and 
may contain a non- virial component to th e line width 
(e.g.- IShen et al.ll2008t iRichards et al.ll2011l ). although it 
is currently unclear if this creates a significant bias in the 
Civ-based mass estimates. We have tried to incorporate 
some systematic errors resulting from a bias in the virial 
mass estimates to the extant that it can be modeled as 
having a simple luminosity-dependence. There were no 
significant differences in our main scientific conclusions 
when we included a simple luminosity-dependent bias. 
In addition, we have tried to mitigate the affects of some 
objects having large systematic errors in the mass esti- 
mates through our use of a Student's i-distribution as a 
model for the measurement errors on the mass estimates, 
which downweights any outliers. 

In general, the number densities of the virial mass es- 
timates obtained before correcting for incompleteness do 
not show discontinuities across redshift bins when switch- 
ing emission lines (Paper I). However, there are a few 
discontinuities in the derived mass and Eddington ratio 
functions when switching between emission lines, sug- 
gesting that small systematics between the virial mass 
estimates of different emission lines may be manifested 
more strongly in our derived BHMF and BHERF. These 
include a discontinuity in the normalization of the BHMF 
when going from H/3 to Mgn (Figure [5]), a discontinuity 
in the BHERF when going from Mgn to Civ (Figure fTO]). 
and differences in the location of the 95 th percentile of 
the distribution of L/LEdd at fixed Mbh inferred from 
the three emission lines (Figure fT2|) . While it may be 
tempting to conclude from this that one or more of the 
emission lines do not give consistent results, this is not 
necessarily the case. As discussed in § I3.2[ the increase in 
number densities in the redshift bins using H/3 is likely at 
least in part caused by a host galaxy contribution to the 
nuclear emission which is not accounted for in the selec- 
tion function, creating an excess of AGN in these redshift 
bins as their nuclear emission gets boosted above the flux 
limit. The discontinuity in the number densities across 
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z ~ 2, corresponding to the shift from Mgn to Civ, also 
corresponds to a transition in the selection function to 
rcdshifts where quasar colors are similar to star colors. 
As discussed above and in § 14.11 this can lead to signif- 
icant systematic error in the incompleteness correction; 
indeed, these redshift bins have the largest statistical er- 
rors as well. 

The only systematic differences between the emission 
lines in the inferred distribution which cannot be ex- 
plained as systematics from incorrect incompleteness cor- 
rection is the anomalous behavior of the 95 th percentile 
of the Eddington ratio distribution at fixed Mbh- This 
quantity is flat with respect to Mbh for the two bins em- 
ploying H/3, increases with increasing Mbh for the red- 
shift bins employing Mgn, and then decreases to a flat or 
approximately flat relationship again for the redshift bins 
employing Civ. Although it has be en argued that Civ is 
less reliable than Mgn or H/3 (e.g. , iBaskin fc Laorll2005t 
iShen et alf2008t iShen fc Liull2012D . in this case it is Mgn 
which gives the discrepant result. It is unclear why this 
is the case, and it may be that this observed trend rep- 
resents real evolution in the joint distribution of Mbh 
and L/LiEdd- However, we consider it likely that uniden- 
tified systematics are at least in part driving this trend. 
Because of this, our current understanding of virial mass 
estimates may not provide enough accuracy for inferring 
how percentiles of the Type 1 quasar joint distribution 
vary in the mass-luminosity plane. 

6. SUMMARY 

In this work we have employed a Bayesian analysis 
method to derive the black hole mass and Eddington 
ratio functions for broad line AGN using a uniformly- 
selected sample from the SDSS DR7. We used more flexi- 
ble models than those in Paper I to test the robustness of 
the key conclusions of Paper I, and found that the main 
results in Paper I remain valid, although the constraints 
are weakened due to the more flexible models used in 
the current work. Our conclusions are summarized as 
follows: 

• The SDSS is < 10% complete at M BH < 3x 1O 8 M 
or L/LiEdd % 0.07, with some variation with red- 
shift. Decreasing the magnitude limits to i ~ 24, 
similar to that of the COSMOS spectroscopic sur- 
veys or the Pan-STARRS medium-deep fields, re- 
duces the mass and Eddington ratio incompleteness 
limits by about an order of magnitude. 

• There are a broad range of Mbh and L/LEdd for 
Type 1 quasars, and there is no evidence for a 
peak in the black hole mass or Eddington ratio 
functions down to the 10% completeness limits of 
the SDSS sample. The number densities of Type 
1 quasars continue to increase toward lower Mbh 
and L/LEdd down to at least Mbh ~ 5 x 10 8 Mq 
and L/LEdd ~ 0.07, respectively. 

• Both the BHMF and BHERF show evidence for 
downsizing. Relative to the peak in the number 
densities at z ~ 2, the number densities of the most 
massive black holes in Type 1 quasars fall off fastest 
toward lower redshift, while the number densities 
of the less massive black holes fall off faster toward 
higher z. This implies that the most massive black 



holes were active first, and shut off their activity 
more rapidly after the peak. 

The number densities of Type 1 quasars are ap- 
proximately constant at z < 2 for L/LEdd ^ 0.05; 
however, the number densities of Type 1 quasars 
radiating at L/Lsdd ^0.1 fall off more rapidly to- 
ward higher redshift at z > 2, possibly reflecting 
a smaller contribution from weaker AGN activity 
toward higher redshift. 

• We constrain the maximum value of L/LEdd in a 
Type 1 quasar to be ~ 3. Therefore, if quasars do 
violate the Eddington limit, they do so only mildly 
and for a short period of time. 

• At low (z < 0.8) and high (z > 2.65) redshifts the 
95 th percentile of the Eddington ratio distribution 
at fixed Mbh is approximately constant with re- 
spect to Mbh- However, at redshift 0.8 < z < 2.65 
the 95 th percentile of p{L/ LEdd\MsH) increases 
toward higher Mbh- This therefore implies that 
at intermediate rcdshifts the shape of the Edding- 
ton ratio distribution changes such that the high 
L/LEdd tail becomes more dominant at higher 
Mbh- At low and high z the shape of the high 
L / LEdd tail of the Eddington ratio distribution is 
independent of Mbh- The redshift dependence of 
this trend is unexpected and may be due to uniden- 
tified systematics among the emission lines used to 
calculate the FWHM-based virial mass estimates. 

• We do not find statistical evidence for the so-called 
"sub-Eddington boundary" in the quasar mass- 
lumino sity plane claimed by ISteinhardt fc Elvisl 
([2010ah . The appearance of such a boundary in 
the "observed" distribution is caused by selection 
effects and errors in the virial BH mass estimates 
(Paper I). This reinforces our early conclusions in 
Paper I that one should not interpret the observed 
distribution directly. 

• Assuming a radiative efficiency of e r = 0.1 and a 
seed black hole mass of Mbh = 10 6 Af Q , the growth 
times for SMBHs in Type 1 quasars having Mbh ^ 
5 x 10 8 M Q are comparable to or longer than the age 
of the universe at z > 1.8. Here, the growth times 
were calculated assuming that SMBHs in a given 
mass bin accrete at a time-averaged rate that is 
equal to the mean Eddington ratio in that mass 
bin. These large growth times imply that prior to 
us observing them as Type 1 quasars, these SMBHs 
experienced a stage of accelerated growth (i.e., with 
higher Eddington ratios). 

• Comparison of the Mbh completeness limits of our 
sample with the local mass function of all SMBHs 
implies that our sample is dominated by SMBHs 
representing the high mass end of the BHMF, 
which reside in what are locally early type galaxies. 
This conclusion in combination with our results on 
SMBH growth times, is consistent with models by 
which SMBHs experience a massive fueling event 
which initiates obscured growth. The black hole's 
growth is self-regulated, persisting until either feed- 
back energy unbinds the obscuring gas or all of the 
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gas is consumed from star-formation, briefly reveal- 
ing the massive black hole as a Type 1 quasar with a 
decaying lightcurve. This same fueling event leaves 
behind a spheroid, placing the SMBH on the Mbh~ 
<7* relationship. Because the SMBHs in our sample 
represent the high mass end of BHMF, this process 
may only be common among this mass range. In 
addition, the long growth times of z < 0.8 Type 1 
quasars with massive BHs and low Eddington ra- 
tios likely represent weaker AGN activity reignited 
by mass loss from evolved stellar populations. 

The combination of our large uniformly-selected sam- 
ple with our powerful Bayesian method represents an im- 
portant contribution to AGN demographic studies, and 
we have used our sample and method in a two-paper se- 
ries to obtain constraints on the optical luminosity func- 
tion, black hole mass function, and black hole Edding- 
ton ratio functions of Type 1 quasars. In addition, we 
have used our sample and method to place constraints 
on the distribution and biases of the virial mass esti- 
mates. The results and methods presented in Paper I 
and in this paper represent important empirical tools 
for understanding black hole growth, for comparison to 
theoretical models, and for planning future surveys. In 
many ways systematic AGN demographic studies with 
respect to Mbh and L/LEdd are just beginning. Further 



improvement will result from using future reverberation 
mapping campaigns to refine our understanding of virial 
mass estimators, as well as applying methods similar to 
our Bayesian technique to current and future deeper sur- 
veys with possibly multiwavelength data. 
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